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ABSTRACT 

The  paper  is  concerned  with  calculating  the  equilib- 
rium configuration  of  a  plasma  confined  by  a  helically  symmetric 
magnetic  field  under  magnetohydrodynamic  assumptions.   This 
entails  solving  an  elliptic  partial  differential  equation  on  a 
domain  with  a  free  boundary.   To  do  so,  we  separate  the  problem 
of  solving  the  partial  differential  equation  from  that  of  com- 
puting the  location  of  the  free  boundary.   A  rectangle  with 
fixed  boundary  is  introduced  as  an  auxiliary  domain  and  the 
differential  equation  is  solved  there.   To  calculate  the  loca- 
tion of  the  free  boundary,  two  variational  problems  are  con- 
sidered.  One  is  the  minimization  of  the  Dirichlet-Douglas 
integral  for  conformal  mapping;  the  other  is  the  minimization  of 
an  integral  connected  with  the  energy  of  the  free  boundary 
problem.   The  method  of  steepest  descent  is  used  to  develop  an 
Iterative  procedure  suggested  by  the  problem  of  minimizing 
these  two  functionals.   Convergence  of  this  algorithm  yields  a 
conformal  mapping  of  the  rectangle  onto  the  original  domain, 
allowing  us  to  express  the  solution  of  the  partial  differential 
equation  there.   The  conformal  mapping  has  the  additional 
property  that  not  only  is  the  solution  of  the  differential  equa- 
tion defined  on  the  original  domain,  but  it  also  satisfies  the 
free  boundary  condition  there. 


The  equations  In  this  one  method  are  approximated  Dy  a 
finite  difference  scheme  which  is  accurate  to  second  order. 
A  program  to  implement  it  has  been  written  for  the  CDC  660O 
computer.   The  program  has  been  run  using  a  variety  of  magnetic 
fields  suggested  by  experimental  machines  which  use  helical 
magnetic  fields  for  plasma  confinement.   The  parameters  have 
been  chosen  so  that  the  ratio  of  plasma  diameter  to  chamber 
diameter  obtained  is  close  to  that  of  the  machine  SCYLLAC. 
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Introduction 

With  the  achievement  of  thermonuclear  fusion  In  the 
hydrogen  bomb,  scientists  began  to  consider  ways  of  controlling 
this  reaction  and  using  it  as  a  power  source.   This  goal 
appeared  to  be  reasonably  easy  to  obtain:   confine  enough 
plasma  at  a  high  enough  density  for  the  necessary  time  and 
fusion  would  occur.   A  magnetic  field  was  chosen  as  the  only 
practical  way  of  containing  the  plasma,  and  several  machines 
based  on  this  principle  were  built.   However,  all  of  them  were 
unsuccessful.   Since  then,  enormous  advances  have  been  made, 
but  even  the  most  successful  device  to  date,  the  Soviet  machine 
TOKAMAK,  is  still  two  or  three  orders  of  magnitude  away  from 
attaining  the  goal  of  a  self-sustaining  thermonuclear  reaction 
(cf.  [2],  [7],  [14]). 

The  question  of  a  mathematical  description  of  the  behavior 
of  a  plasma  in  a  magnetic  field  has  not  been  neglected.   In 
this  paper,  we  study  one  aspect  of  the  problem.   We  adopt  a 
magnetohydrodynamic  model  to  describe  the  behavior  of  a  plasma 
and  we  specify  a  numerical  method  for  the  computation  of 
figures  of  equilibrium  when  the  magnetic  field  is  helically 
symmetric.   This  configuration  is  applicable  to  the  machine 
SCYLLA  IV  and  qualitatively  applicable  to  future  toroidal 
machines,  such  as  SCYLLAC. 


In  addition  to  its  physical  applications,  this  problem 
is  also  interesting  mathematically.   It  requires  the  solution 
of  an  elliptic  partial  differential  equation  on  a  domain  with 
a  free  boundary,  a  problem  which  has  long  been  of  interest  in 
the  field  of  partial  differential  equations  (cf,  [13]). 

Previous  research  into  the  numerical  solution  of  free 
boundary  problems  has  been  done  by  E.  Bloch  [5,6],  who  con- 
siders partial  differential  equations  of  the  form 


^xx   ^yy 


where  the  dots  stand  for  first  order  terms.   He  solves  these 
problems  by  introducing  a  fixed  auxiliary  domain  in  the  shape 
of  a  rectangle  and  then  using  the  method  of  steepest  descent  to 
obtain  a  conformal  mapping  from  the  auxiliary  domain  to  the 
physical  plane.   After  that,  it  is  only  necessary  to  solve  the 
appropriate  Dirichlet  problem  on  the  rectangle. 

In  this  paper  we  extend  Bloch 's  work  and  study  more 
general  elliptic  partial  differential  equations 


ail/      +  2bii/       +  cTp       +  All/    +  etl/    +  ff  =   g,    , 
^xx  ^xy        ^yy        ^x        ^y        ^        ^    ' 


where  a,  b,  c,  d,  e,  f  and  g  are  functions  depending  on  the 
independent  variables  x  and  y.   By  allowing  the  second  order 
terms  to  have  variable  coefficients  we  increase  the  complexity 
of  the  problem.   First,  the  mapping  we  obtain  is  conformal  with 
respect  to  the  metric 


2       2  2 

ds  =  ady  -  2bdxdy  +  cdx 


Instead  of  the  usual  Euclidean  norm.  Second,  the  physical 
domain  where  we  solve  the  problem  is  doubly-connected.  This 
fact  implies  that  our  auxiliary  domain  cannot  be  a  rectangle, 
but  must  be  a  ring.  Furthermore,  we  allow  the  fixed  boundary 
to  be  of  arbitrary  shape.  In  our  case  the  algorithm  does  not 
have  a  unique  solution  and  it  is  necessary  to  devise  a  scheme 
to  correct  this  deficiency. 

The  description  of  the  computation  will  be  accomplished 
in  four  stages.   In  the  first  chapter  the  physical  problem  is 
defined,  while  in  the  second  an  equivalent  variational  formula- 
tion is  given.   In  the  third  chapter  we  discuss  a  method  of 
obtaining  a  conformal  mapping  by  variational  methods.   The 
fourth  chapter  combines  these  two  variational  problems  to  derive 
the  actual  algorithm  used  in  performing  the  calculations.   In 
the  final  chapter  the  results  of  some  calculations  are  given 
and  their  accuracy  is  discussed.   The  actual  computer  program 
used  is  included  as  Appendix  C. 


Chapter  1  -  The  Physical  Problem 

1. 1   Introduction 

Magna tohydrodynamlcs  Is  one  of  the  models  used  In  plasma 
physics  for  describing  the  behavior  of  an  Ionized  gas,  or 
plasma,  contained  by  a  magnetic  field.   It  Is  quite  useful, 
since  It  combines  reasonable  physical  accuracy  with  mathematical 
simplicity.   In  the  magnetohydrodynamlc  approach,  the  plasma  Is 
considered  to  be  a  classical  macroscopic  conducting  fluid, 
characterized  by  its  density  p,  its  entropy  S  and  its  flow 
velocity  u.   The  pressure  p  is  then  a  function  of  p  and  S.   The 
magnetic  field  B  and  the  current  J  which  Interact  with  this 
plasma  exert  a  Lorentz  force  J  x B  on  the  plasma.   A  set  of  non- 
dissipatlve  fluid  equations  coupled  with  Maxwell's  equations 
through  the  Lorentz  force  term  and  Ohm's  law  are  then  used  to 
obtain  the  necessary  formulas  to  describe  the  behavior  of  the 
plasma. 

We  make  the  further  assumption  that  the  plasma  is  perfectly 
conducting.   Also,  the  charge  density  and  the  displacement 
current  are  neglected.   In  this  case.  Ohm's  law  for  a  perfect 
conductor  states  that 

(1)  E  =  B  xu  , 

where  E  is  the  electric  field  vector.   If  we  adopt  a  frame  of 
reference  moving  with  the  fluid,  then  E  must  vanish.   Since  we 
are  ignoring  the  displacement  current,  we  can  express  J  in  terms 


of   the   magnetic   field   by 

J    =  VxB    ; 

it  is  now  possible  to  eliminate  all  of  the  electromagnetic 
variables  except  B  from  Maxwell ' s  equations.   The  integral  forms 
of  these  equations  are  then 


(2)  /  B-vda  =  0 
and 

(3)  ^  J  B-vda  +p(Bxu)d^  =  0  , 

S      aD 


where  D  is  an  arbitrary  domain  with  inner  normal  v  and  surface 
element  da  and  S  is  an  open  surface  contained  in  D.   Since  the 
frame  of  reference  is  moving  with  the  fluid,  the  domain  D  can 
vary  with  time.   The  derivative  -rr-   which  occurs  in  equation  (3) 
is  the  Lagrangian  time  derivative  defined  by 

dt     at 

Equation  (3)  is  just  Faraday's  law  in  integral  form.   Its 
validity  allows  us  to  interpret  equation  (2)  as  an  initial  condi- 
tion, instead  of  an  equation  of  motion. 

Turning  now  to  fluid  mechanics,  we  set  down  three  equations 
for  the  conservation  of  mass,  momentum  and  energy.   The  equation 
for  the  conservation  of  mass  is 


/  ^  '^^  ^  f  ^p^^"^^  =  °  • 

D  SD 


The  divergence  theorem  can  be  applied  to  the  surface  Integral 
to  bring  it  into  the  more  convenient  form 

(h)  f  i^  +   div  (pu)}dV  =  0  . 


)t 
D 


The  equation  for  the  conservation  of  momentum  is  derived 
from  Newton's  second  law  of  motion.   The  two  forces  acting  are 
the  mechanical  force  -Vp  and  the  Lorentz  force  Vx  (VxB). 
This  equation  therefore  has  the  integral  form 

(5)         ^  j  pudV  +  j  VpdV  -  j  (JxB)dV  =  0  , 

D        D        D 

where  -s-p  is  again  the  Lagranglan  time  derivative- 

The  third  conservation  law  derived  from  fluid  mechanics 
is  for  the  entropy  S.   It  is  assumed  that  the  flow  is  adlabatic 
and  that  the  fluid  acts  like  a  perfect  gas.   Then,  under  these 
conditions 


'6'  al=° 


The  system  of  equations  (5),  (4),  (5)  and  (6),  together 
with  the  initial  condition  (2),  is  known  as  the  Lundquist  equa- 
tions.  These  equations  can  also  be  expressed  in  the  more  usual 
differential  form 


§1+  (V.u)B-  (B.V)u  :=  0  , 


P  §^  +  Vp  +  Bx  (VxB)  =  0  , 


p  dt      ^    u  , 


and 

dS 


dt  -  0  ' 


on  the  domain  D,  along  with  the  Initial  condition 

V.B  -  0  . 


1. 2  Derivation  of  the  Equilibrium  Equations 

In  this  paper,  we  shall  only  be  concerned  with  determining 
figures  of  static  equilibrium  In  magnetohydrodynamlcs  for  a 
perfectly  conducting  fluid  which  Is  contained  by  a  magnetic 
field.   In  this  case  there  Is  no  flow,  and  therefore  the  velocity 
must  vanish.   Several  consequences  follow  from  this  assumption. 
First,  as  a  result  of  Ohm's  law  (l),  the  electric  field  E  must 
be  zero.   Next,  equation  (4)  states  that 

(7)  ^  =  0 
^' ^  dt   ^  ' 

and  so  p  must  be  a  constant.   The  Lundquist  equations  thus 
reduce  to  the  simpler  system  of  equations 

(8)  JjB-vdo  =   0 

^D 
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and 


(9)  if  ^P  +  ^)^-   (B-v)B]da  =  0 


In  these  conditions  for  equilibrium,  equation  (8)  is  no  longer 
to  be  treated  as  an  initial  condition. 

These  integral  relations  are  basic  to  the  formulation  of 
the  problem.   Under  assumptions  of  smoothness  they  will  lead 
to  the  necessary  differential  equations,  otherwise  to  the 
appropriate  jump  conditions.   We  first  consider  the  smooth  case 
where  the  domain  D  lies  entirely  in  one  medium  and  both  B  and  p 
are  continuous.   The  integral  formulation  then  yields  the 
differential  equations 

(10)  V.B  =  0 
and 

(11)  Vp  +  Bx(VxB)=:0. 

On  the  other  hand,  D  can  lie  across  an  interface  r  which 
separates  two  distinct  mediums.   Both  integral  equations  (8) 
and  (9)  remain  valid  in  this  case,  even  if  the  domain  D  is  taken 
arbitrarily  small.   Under  the  assumption  that  B  is  bounded,  we 
obtain  the  fact  that  both  integrands  must  be  continuous  across 
the  discontinuity  surface.   Hence,  both  B-v  and  -p- +  p  must  be 
continuous  across  r. 
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1. 5  Description  of  the  Physical  Problem 

We  shall  now  describe  the  equilibrium  problem  which  will 
be  considered  In  this  paper.   Let  P  be  a  plasma  which  Is 
surrounded  by  a  vacuum  D.   The  common  boundary  separating  the 
plasma  and  the  vacuum  will  be  denoted  by  T,    and  the  outer 
boundary  of  D  by  C.   We  shall  assume  further  that  the  location 
of  C  Is  fixed,  while  that  of  r  Is  defined  by  the  equilibrium 
position  of  the  plasma. 

Inside  P,  neither  current  nor  magnetic  field  Is  present. 

(P) 
Hence  B^  ' ,    the  magnetic  field  in  the  plasma,  is  zero.   Under 

these  conditions,  the  differential  equations  (10)  and  (11)  yield 

the  result  that  the  pressure  p  reduces  to  a  constant.   Inside  D, 

the  pressure  is,  of  course,  zero.   Hence  B^   ,  the  magnetic 

field  in  the  vacuum,  satisfies 

(12)  V.B^^^  -  0 
and 

(13)  VxB^^^  -  0  , 

the  appropriate  versions  of  equations  (8)  and  (9)-   Moreover, 
a  consequence  of  equation  (13)  is  that  the  current  J^  '  in  the 
vacuum  must  be  zero. 

At  the  common  boundary  r  between  P  and  D,  the  jump  condi- 
tions which  we  derived  in  the  last  section  yield  the  boundary 
conditions 


(14)  B^^^v  =  0 


and 

(15)  I  [b'^'i^  =  P 


Similar  reasoning  at  the  outer  boundary  C  of  the  domain  D  yields 
the  simpler  boundary  condition 

(16)  B^^^V  =  0  . 

We  can  derive  the  free  boundary  condition  (15)  in  another 
way.   Even  though  no  current  exists  in  either  the  plasma  or  the 
vacuum,  the  discontinuity  in  the  magnetic  field  at  the  interface 
r  can  cause  a  current  to  flow  there.   This  surface  current  K  is 
given  by 

K  =   VXB^^)  , 

b(v) 

while  the  average  value  of  B  on  r  is  — ^ — .   Hence, the  Lorentz 

rB(V)n2 
force  at  the  interface  r  is  — p  •'  v,  and  we  may  interpret 

equation  (15)  as  stating  that  for  a  figure  of  equilibrium  to 
exist,  the  fluid  pressure  of  the  plasma  must  balance  the  mag- 
netic pressure  exerted  by  the  magnetic  field. 

f  Pi 

Since  B^  ''  is  identically  zero,  we  may  drop  the  super- 
scripts from  our  notation  and  refer  unambiguously  to  the  magnetic 
field  B^"^^  defined  on  D  as  B. 

In  addition  to  the  differential  equations  (12)  and  (13) 
and  the  boundary  conditions  (14),  (15)  and  (l6),  the  complete 
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free  boundary  problem  requires  the  specification  of  an  appro- 
priate number  of  periods.   The  exact  number  is  determined  by 
the  geometry  of  the  problem  (cf.  [4]). 

These  periods  may  be  specified  in  two  ways,  either  as 
fluxes  or  as  currents.   A  flux  0„   is  defined  by 


(17)  03  =ffB'<io   , 


S 
where  S  is  a  surface  whose  boundary  lies  in  the  boundary  of  D, 
while  a  current  In  is  defined  by 


(18)  I_g  =(7)B-ds  , 

where  i  is  a  closed  loop. 


1. 4   Comments  on  the  Problem 

Of  even  greater  importance  than  the  computation  of  figures 
of  magnetohydrodynamic  equilibrium  is  the  determination  of  their 
stability.   While  we  will  not  discuss  the  stability  of  the 
solutions  obtained,  the  method  employed  will  give  some  indica- 
tion of  stability.   The  numerical  scheme  is  based  on  the  method 
of  steepest  descent  and  on  the  idea  of  making  the  energy 

(19)  E=liF«'^^ 

D 

a  minimum.   The  method  of  steepest  descent  can  be  pictured  as 
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following  the  steepest  path  downward  until  a  minimum  is 
attained.   If  the  energy  is  not  a  true  minimum,  this  will 
disturb  the  workings  of  the  numerical  method  and  prevent 
convergence . 

From  a  mathematical  point  of  view,  the  free  boundary 
problem  is  quite  difficult  to  solve.   The  problem  is  both 
three-dimensional  and  contains  a  free  boundary.   However,  if 
we  assume  that  B  has  a  symmetry,  then  the  number  of  independent 
variables  may  be  reduced  to  two  and  the  whole  problem  becomes 
less  complex.   Moreover,  we  can  use  the  condition  V-B  =  0  to 
Introduce  a  flux  function  and  simplify  the  problem  even  further. 

In  this  paper  we  shall  make  the  assumption  that  the  mag- 
netic field  B  is  helically  symmetric.   This  assumption  allows 
us  to  make  the  mathematical  simplifications  described  in  the 
last  paragraph  and  also  corresponds  to  physical  reality. 
Certain  machines  for  controlled  thermonuclear  power,  such  as 
the  Stellerator,  have  been  constructed  in  the  shape  of  a  long 
thin  torus  with  helical  windings  wrapped  around.   These  windings 
give  rise  to  a  helically  symmetric  magnetic  field. 

In  the  next  sections  we  will  give  a  precise  definition  of 
helical  symmetry  and  determine  what  conditions  a  magnetic  field 
must  fulfill  in  order  to  be  both  helically  symmetric  and 
satisfy  the  free  boundary  problem. 
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1-5  Helical  Symmetry 

To  define  helical  symmetry,  we  start  with  the  framework 
of  the  cylindrical  coordinate  system  (r,9,z).  In  this  system 
the  vector  B  is  written  as 


(20)  B  =  B^e^+B^eg+B^e^  , 


where  e  ,  e^  and  e   are  the  canonical  basis  vectors  and  B  ,  Bg 
and  B   the  components  in  these  directions.   We  will  say  that 
the  vector  B  is  helically  symmetric  if  B  depends  on  only  the 
two  independent  variables  r  and  ^,    where  r  is  the  same  variable 
as  in  equation  (20)  and 

(21)  (f)  =  e  -  kz  . 

In  this  paper  we  shall  be  interested  only  in  those  solutions  of 

the  free  boundary  problem  which  have  this  property,  and  the 

present  section  will  be  devoted  to  formulating  the  conditions 

such  a  helically  symmetric  vector  must  satisfy  in  order  to  be 

a  solution. 

Since  B  depends  on  only  two  variables,  and  not  three,  it 

will  be  convenient  to  introduce  a  new  set  of  basis  vectors 

which  distinguishes  the  direction  in  which  B  remains  constant. 

We  therefore  construct  a  new  set  of  basis  vectors  e  ,  ei  and  e^.. 

r'   (J)      C 

The  first  two  vectors  will  span  the  subspace  where  B  varies, 
while  Q^   will  lie  in  the  ignorable  direction  where  B  does  not 
vary.   The  vector  e   is  the  same  basis  vector  as  in  cylindrical 
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J 


coordinates.   The  vector  ei  is  defined  in  terms  of  the  original 
set  of  basis  vectors  as 

ep.  -  kre^ 

^  a 

where 

(23)  a  =  1  +k^r   . 


Equation  (22)  defines  the  second  direction  in  which  a  helically 
symmetric  function  can  vary.   To  see  this,  let  ^(r,0,z)  be  a 
helically  symmetric  function.   The  gradient  of  ip   expressed  in 
cylindrical  coordinates  is 

\      J  ^   ^r  r    r   9   z  z 

However  ip   is  helically  symmetric  and  depends  only  on  the  variable 

<t)  =  0  -  kz  . 
Hence  we  have 

(25) 
and 


^f   _  hip 


(26)  M  .  _k|?  . 

dz       d(p 

Inserting  these  relations  into  equation  (24),  we  see  that 


(27)  VV/  =  ^^e^  +  T"  ^^ 


2    f  Qd  -  kre 

a   ,  /  9     z 


2 
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Therefore  the  two  vectors  e^  and  ei  span  the  space  of  Vtl/. 

The  vector  e^  then  lies  in  the  perpendicular  direction 

In  which  ^  remains  constant.   We  use  standard  methods  of  linear 

algebra  to  construct  this  vector  and  obtain 


(28)  e,  -     ^  '^ 


kreg  +  e_ 


Thus  we  have  a  set  of  orthogonal  basis  vectors  with  which 
to  discuss  helical  symmetry.   The  vector  B  is  written  in  terms 
of  these  new  basis  vectors  as 

B  =  B  e  +Biei  +B^e^    , 
r  r   9  9   C  C 

where  the  components  Bi  and  B^   are  defined  by 

B^  =  Bg  -  krB^ 

and 

B^  =  krBo  +  B   . 

Q     y   z 

Once  these  transformations  are  defined,  we  can  write  all 
the  usual  vector  operators  in  terms  of  these  new  basis  vectors. 
A  complete  list  will  be  found  in  Appendix  A.   The  free  boundary 
problem  now  consists  of  solving  the  equations 

.  S(rB  )    -,  SB  I 
and 
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=   0 


on  the  interior  of  D,  with  the  fixed  boundary  condition 


(51) 


B'V  =  0 


and  the  free  boundary  condition 


(32) 


2    2 

a    a 


-  2p 


1. 6  The  Flux  Function 

Since  the  magnetic  field  B  now  depends  on  only  two  inde- 
pendent variables,  it  is  possible  to  introduce  a  flux  function 
Tj/   and  reformulate  the  problem  In  terms  of  this  function.   We 
take  ^(r,(}))  to  be  a  continuous,  helically  symmetric  function 
defined  on  the  domain  D  by 


(33) 

and 

(34) 


rB 


-B 


^  • 


This  definition  of  tl/   clearly  satisfies  equation  (29),  and  we  can 
express  the  vector  B  in  terms  of  ?//  as 


B 


=  B^e^  -  e^  X  Vi// 
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Equation  (30)  then  leads,  on  one  hand,  to  the  elliptic  partial 
differential  equation 

and,  on  the  other,  to  the  condition 

B^  =  constant  . 

The  boundary  conditions  (3I)  and  (32)  are  also  expressed  easily 
in  terms  of  ip.      The  boundary  condition 

B-v  =  0 

Is  equivalent  to  the  condition 

■>p  =   constant  , 

while  the  free  boundary  condition  (32)  becomes 

\  r     o  o     ■' 

It  can  be  shown  that 

B-V^  =  0  ; 

the  magnetic  field  is  perpendicular  to  the  gradient  of  ^.   Hence 
the  level  lines  of  f   are  parallel  to  the  magnetic  field,  and  we 
verify  that  t//  is  a  flux  function  for  B. 
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The  appropriate  periods  for  this  formulation  of  the 
problem  are  given  in  terms  of  the  fluxes  (cf.  equation  (17)) 
Two  periods  must  be  specified.   Possible  periods  are 


0^   =  B^rdrde 


D 

the  flux  in  the  z  direction, 


^o-Ij 


Bgdrdz 
D 

the  flux  in  the  0  direction,  and  H  the  difference  in  the  value 
of  ■!//  between  the  fixed  and  free  boundaries.  These  three  terms 
may  also  be  expressed  in  terms  of  the  flux  function  ?//  as 


2 

(37)  ^z=jj( — v^;^'-<^<'  • 

D       "^ 


(38)  00  -  [[(^   -  ^  /drd^ 

^P^  ">  ka    a   / 

and 


(39)  H  =  ^^  -^P  , 


where  ^^  is  the  constant  value  of  tp   on  the  fixed  boundary  C  and 
^P  is  the  constant  value  on  the  free  boundary  r.   Since  the  flux 
function  i^  is  defined  up  to  a  constant,  we  normalize  ^  so  that 


^P  =  0  . 
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All  flux  functions  which  we  will  consider  in  this  paper  will 
be  assumed  to  have  this  property.   Thus,  relation  (39)  is 
equivalent  to  the  boundary  conditions 

^  =  0 

on  r  and 

^  =  H 

on  C. 

The  three  periods  0  ,  0g  and  H  are  related  by  the  identity 

(^0)  H  =  4  (00  +0^)  • 

Therefore,  if  we  specify  any  two  of  the  set  0  ,  0^  and  H,  equa- 
tion (40)  allows  us  to  determine  the  value  of  the  third  one.   In 
addition,  the  constant  B^,    whose  value  has  not  yet  been  given, 
is  now  defined  implicitly  by  either  equation  (37)  or  equation  (38) 

The  complete  free  boundary  problem  can  now  be  expressed 
in  terms  of  f.      On  the  domain  D  we  wish  to  find  a  function  ipir,^) 
which  satisfies  the  partial  differential  equation  (35)  on  the 
interior  of  D,  the  free  boundary  condition  (36)  on  the  boundary 
r  of  D  and  the  three  relations  (37),  (38)  and  (39).   As  we  said 
above,  only  two  of  the  constants  0  ,  0q   and  H  must  be  specified; 
the  third  is  determined  by  equation  (40).   This  formulation  of 
the  problem  is  the  one  we  will  solve  numerically. 


19 


1.7  The  Potential  Function 

While  the  expression  of  B  in  terms  of  a  flux  function  is 
the  one  we  will  use  when  solving  the  free  boundary  problem,  it 
is  also  of  interest  to  have  an  expression  for  B  in  terms  of 
the  potential  function  O.   In  this  case  we  start  from  equation 
(30)j  the  condition  that  V  xB  =  0.   Written  in  cylindrical 
coordinates,  this  is 

/  .  hB  SBnN      /  SB    hB 

(41)     VxB=e   1^  -  ^;  +  e.  U-^  -  -^j 
^   '  r  \  r  dS    Sz  y    9  \dz  dv     j 

Clearly,  each  component  of  this  vector  must  be  zero.   The  first 
component  of  the  vector  states  that 

4  (B,+krBg)  =  0  , 

while  a  linear  combination  of  the  last  two  components  shows  that 


^ 


(B  +krBn)  =  0 


Therefore,  the  component 


B^  =  krBg  +  B^ 


must  be  a  constant.   We  now  define  the  continuous,  helically 
symmetric  function  C){t,^)    on  the  domain  D  by 

(^2)  1^  =  B^ 

^   '  dr    r 


eo 


and 


,,  ,  So 

2 

rB  1  +  kr  B,, 

=      g     . 

o 

If  we  now  define 

then  clearly  O  satisfies  equation  {30),    and 

(44)  B  =  VO(r,0-  kz)  +B^z  . 

Note  that  the  potential  function  for  B  is  not  helically  symmetric, 
even  though  the  cylindrical  components  of  B  are. 

Equation  (29)  leads  now  to  the  elliptic  partial  differential 
equation 

(45)  ^  (rO^)^  +^"H  =  °  • 

In  addition,  the  boundary  condition 

B-v  =  0 
implies  that 

2 

(46)  rO  i^  -  —  Oir.  +krB;.r   =  0 

on  both  the  free  boundary  r  and  the  fixed  boundary  C. 

The  appropriate  periods  for  this  formulation  of  the 
problem  are  currents  (cf.  equation  (l8)).   The  three  possible 
periods  are 
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^z  =,^  (B^dr+rBgd0) 


the  current  in  the  (r,0)  plane, 


Ip   -IJ   (B^dr  +B^dJ 


the  current  in  the  (r,z)  plane  and  B^,    the  component  of  the 
magnetic  field  B  in  the  ignorable  direction.   As  in  the  case 
for  the  flux  function,  there  exists  an  algebraic  relation 
between  these  three  periods  given  by 


(47)  B^=lF(^z-Ie) 


Therefore  it  is  sufficient  to  specify  any  two  of  I  ,  Ig  and  B^ . 
The  functions  f   and  O  are  conjugate,  in  a  sense,  and  the 
relations  between  them  are  given  by 


(48)  ^  =  _  o_  Oi  +  krB» 
^      ^  ^r     r   (f)     (^ 

and 

(49)  f^   -rO^  . 


1.8  The  Scaling  of  B^  and  k 

The  parameters  B^  and  k  may  be  scaled  out  of  the  partial 
differential  equation  (35).   If  we  first  make  the  substitution 

T  =  kr 
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then  our  equation  becomes 


k^Tp 


^^7? 


=  0  . 


(1+  T  ) 


We  now  Introduce  the  function  ^{t,(|))  defined  by 


i>{^M   =  4-  t(T,(t))  . 


B 


This  substitution  further  simplifies  the  equation  to 


(50) 


^\ 


+ 


1  /' 


T^ 


27? 


"^V  1+  T  h  (1  +T^) 


=  0 


The  boundary  conditions  on  i\i   can  be  easily  derived  from  those 
on  ill.      On  the  fixed  boundary  C  we  have 

^        kH 


^  = 


B  i» 


while   on   the   free   boundary  r   we  have 


and 


ii  =  Q 


B 


^ 


li 


•72 


+ 


+ 


1+   T  1  +  T 


+    2p     . 


1. 9  An  Exact  Solution  of  the  Problem 

While  the  general  solution  to  this  problem  in  closed  form 
is  not  known,  we  do  know  the  solution  for  the  case  when  the 
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fixed  boundary  C  Is  a  circle  whose  center  is  at  the  origin  of 
the  coordinate  system.   To  see  this  more  easily,  let  us  examine 
the  problem  when  it  is  expressed  in  terms  of  cylindrical  co- 
ordinates.  In  this  case  there  are  two  simple  solutions.   The 
first  is  a  magnetic  field  lying  parallel  to  the  z-axis. 


B  =  ae   , 


where  a  is  an  arbitrary  constant.   This  solution  leads  to  the 
flux  function 


(51)  ^(r)  =  I  B^kr^  , 


where 


B,.  =  a  . 


The  second  solution  in  cylindrical  coordinates  is  the 
logarithmic  field 

B  =  ^  eg  , 

where  p  is  an  arbitrary  positive  constant.   Elementary  algebra 
shows  that  this  is  equivalent  to  the  flux  function 


B  lo 
(52)  ^(r)  =  -  ^  log  r  , 


k 


where 


B^  =  Pk  . 
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since  equation  (55)  Is  inhomogeneous,  we  cannot  add 

equations  (51)  and  (52)  to  obtain  a  general  solution.   However, 

we  can  combine  them  In  a  weighted  average.  Therefore,  the 
general  solution  for  this  case  Is 

B^kr^         B» 

(53)      ^(r,({))  -  7  -^^^—  +  (7-1)  ^  log  r  +  5  , 

where  7  and  5  are  arbitrary  constants.   Equation  (53)  will  be 

used  to  check  the  accuracy  of  the  numerical  solutions  which  we 
will  obtain. 
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Chapter  2  -  A  Variational  Approach 

2. 1   Introduction 

In  the  first  chapter  we  have  formulated  the  complete 
physical  problem  which  we  will  solve.   However,  in  order  to 
derive  the  algorithm  which  we  will  use  for  the  problem,  it  will 
be  helpful  to  introduce  an  associated  variational  problem. 
Once  this  problem  has  been  defined,  we  will  compute  the  first 
variation  and  show  that  it  leads  to  the  free  boundary  condition 
(1.15). 

To  this  end,  in  the  next  section  we  will  prove  Dirichlet's 
principle  for  the  flux  function  ^.  Furthermore,  in  Section  2,3  we 
will  introduce  an  appropriate  comparison  function  and  use  the 
Dirichlet's  principle  we  have  just  derived  to  obtain  an  upper 
bound  for  the  first  variation.   An  equivalent  variational 
problem  expressed  in  terms  of  the  normalized  potential  function 
will  be  studied  in  the  last  section.   Here  a  similar  procedure 
will  be  used  to  derive  a  lower  bound  which  agrees  with  the  upper 
bound  to  first  order.   We  will  prove  Thomson's  theorem  for  the 
normalized  potential  function  and,  by  choosing  a  new  comparison 
function,  compute  the  lower  bound  and,  hence,  the  first  varia- 
tion. 

The  procedure  we  will  follow  has  other  applications  in 
addition  to  computing  the  first  variation.   It  can  also  be  used 
to  analyze  the  second  variation  and  thus  to  determine  the 
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stability  of  the  solution,  a  question  of  great  practical  impor- 
tance.  In  fact,  the  same  methods  that  we  will  employ  have  been 
used  (cf.  [3,12])  to  show  the  stability  of  the  cusped  geometry 
model  in  magnetohydrodynamics. 


2. 2  A  Dirichlet  Principle  for  the  Flux  Function 

We  can  formulate  Dirichlet 's  principle  for  the  flux  func- 
tion -ip.      To  recapitulate  briefly,  the  solution  of  the  free 
boundary  problem  is  a  helically  symmetric  vector  B  which  satis- 
fies equations  (1.29)  and  (I.30)  on  the  interior  of  D,  equa- 
tions (1.31)  and  (1.32)  on  the  boundary  of  D  and  has  two  periods 
specified.   Moreover,  this  magnetic  field  B  has  associated  with 
it  an  energy 


(1)  E  =  "I  //  B^rdrd(})  . 


D 

As  was  shown  in  Section  1.6,  we  can  introduce  a  flux  function  f 
for  the  magnetic  field  B  that  satisfies  the  partial  differential 
equation  (I.35),  the  free  boundary  condition  (I.36)  and  the 
fixed  boundary  conditions 

^  =  H 
on  C  and 

^  =  0 

on  r.   The  two  periods  mentioned  above  can  be  any  two  of  the 
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set  (1.57),  (1.38)  and  (l.39)«   However,  by  means  of  relation 
(1.40)  It  Is  possible  to  determine  the  value  of  the  third  one. 
Therefore,  we  know  the  value  of  the  constant  H  and  the  flux 

(2)  0^  =JJ  B^rdrd^    . 

D 

For  purposes  of  comparison,  let  us  Introduce  on  the  same 
domain  D  another  magnetic  field  denoted  by  B  .   This  new  field 
Is  similar  to  B  In  that  It  satisfies  equations  (I.29)  and  (I.5I) 
and  differs  from  B  In  that  It  Is  not  required  to  satisfy  equa- 
tions (1.30)  and  (1.32).   In  addition,  the  energy  associated 
with  B*  Is  defined  as 

(3)  E*  =  I^B^^rdrdc})  . 

D 

Proceeding  as  we  did  for  B,  we  can  use  equation  (I.29)  to 
Introduce  a  flux  function  Tp*  for   B^.   Unlike  ^,  f*   does  not 
satisfy  equation  (1-35)  since  equation  (I.30)  does  not  hold 
for  B*.   On  the  other  hand,  the  fact  that  B  and  B*  have  the 
same  periods  Implies  that  tp*   must  have  the  same  fixed  boundary 
conditions  as  f, 

^*  =  H 


on  C  and 


Tp*  =  0 


on  r.   Also,  the  flux  0     of  ip*   Is  the  same  as  that  of  ip. 
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Dlrichlet's  principle  then  states  that 


(4)  E  <  E 


* 


,* 


for  all  such  magnetic  fields  B  .   We  perform  the  proof  in  terms 
of  the  two  flux  functions  -^  and  f* .      First,  we  express  the  two 
energies  in  terms  of  their  respective  flux  functions,  obtaining 

(5)  E  =  |j7(|.!j.4)-^-''* 

and 

1  rrf'^h   ,  ^J  ,  ^_^ 

r     a     a' 


(6)  E*  =^     -1-  +  -^  +  -^  ^rdrd^ 


2  J  J  \  ^T       ^^  2 

D 


Also,  the  flux  0  may  be  expressed  in  terms  of  either  of  the 
two  flux  functions,  as 


V 

or  as 


OP  kr^  +  B»'n 


D   ■♦    ^ 


kr^^  +  B* 


^.  -  JJ  { -V^  1  -^-^^  ■ 

We  now  define  the  two  functions  f[v,^)    and  g(r,i{))  by 

f(r,(t))  =  ^""(r,.!))  -  ^(r,(t)) 

and 

g(r,(t))  -  B^(r,(l))  -  B^  . 

Since  ?//  and  ip*   both  have  the  same  periods,  it  follows  that 


tj/  =   f* 
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on  the  boundaries  C  and  r.   Hence, 


f(r,<t))  =  0 


there.   Expressing  the  energy  E   in  terms  of  ip ,    B^,  f  and  g, 
we  see  that 


(7) 


E   =  E  + 


D 


^  ^  r^/r  ^  ^^^^ 


+ 


o 


drdij) 


-i 


D 


-^  +  -2  + 
r    a 


a 


rdrdcf) 


Furthermore,  the  condition  that  the  flux  j6      is  identical  for  B 
and  B*"  imposes  the  requirement  that 


D 


/krf  +gN 

f ^ —  I  rdrdc})  =  0 

\    a 


If  we  include  this  equality  in  equation  (7)  and  neglect  the 

2 


positive  term 


D 


2     2 

ff    f 


^ 


+  ^  +  s 


rdrd(J3,  we  obtain 


(8) 


E   >■  E  + 


D 


V^ 


+ 


VTp   f    B.kr^f 
^r  r    c    r 


a 


drd(}) 


By  applying  the  divergence  theorem  to  the  term  on  the  right  and 
by  making  use  of  the  facts  that  ip   satisfies  equation  (1.35)  on 
the  interior  of  D  and  that  f  =  0  on  the  boundary  of  D,  we  get  the 
desired  result, 

E  <  E*  . 
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2. 3  The  Derivation  of  an  Upper  Bound  for  the  First  Variation 

We  will  now  use  variational  methods  to  derive  the  free 
boundary  condition.   In  the  previous  section  we  showed  that  on 
a  given  domain  D  the  solution  could  be  characterized  as  an 
extremal  within  an  appropriate  family  of  magnetic  field  vectors, 
but  we  said  nothing  about  the  shape  of  D.   Similar  methods  will 
now  be  used  to  define  this  shape  and  to  derive  the  free  boundary 
condition 

(9)  ^  =  P 

which  must  be  satisfied  at  each  point  of  r .   We  will  define  a 
family  of  admissible  domains,  each  having  a  unique  magnetic 
field  defined  on  it.   The  domain  D  on  which  the  free  boundary 
problem  is  solved  will  then  be  characterized  as  the  one  which 
minimizes  the  Dirichlet  integral  defining  the  magnetic  energy. 
To  be  more  precise,  we  will  define  a  variation  of  D  onto 
the  domain  D"^.   This  is  done  by  introducing  the  notion  of  an 
interior  variation  of  D.   Let  S-,  (r,c|))  and  Sp(r,(t))  be  a  pair  of 
functions  possessing  plecewise  continuous  partial  derivatives 
of  a  sufficiently  high  order  in  some  neighborhood  of  D.   For 
sufficiently  small  e,  the  transformation 

(10)  r*  =  r+  eS^(r,(i)) 

(11)  h*   =  <t>+  eS^(r,<t)) 
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maps  D  onto  the  nearby  domain  D*  in  a  one-to-one  fashion.   The 
mapping  defined  by  equations  (10)  and  (11)  is  known  as  an 
interior  variation  since  it  is  defined  on  the  interior  as  well 
as  the  boundary  of  D  (cf.  [12]).   We  place  only  two  constraints 
on  this  variation.   First,  the  area  of  D*  must  be  the  same  as 
that  of  D.   Therefore 

(12)  //  v*(ir*6^*  =    //  rdrdif)  . 

D*  D 

Second,  the  location  of  the  fixed  boundary  C  of  the  domain  D 
must  be  identical  with  that  of  the  fixed  boundary  C*  of  D*. 
Hence, 

(13)  S^(r,<t))  =  S2(r,(t))  =  0 

for  all  points  {r,\))    lying  on  the  boundary  C. 

The  free  boundary  r*,  however,  is  not  in  the  same  location 
as  r.   It  is  an  infinitesimal  distance  6v  away  from  r,  where  5v 
is  measured  along  the  inner  normal  from  r.   It  is  understood 
that  if 

5v  <  0 

then  the  displacement  is  in  the  direction  of  the  outer  normal. 
Moreover,  we  can  express  6v  in  terms  of  the  interior  variations 
Sn  and  Sp  by 


(1^)         6v  =  E 


S,(r,^)  If  .  r%(r,<>)  |i 
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Continuing  with  our  construction,  we  Introduce  on  each 
domain  a  unique  magnetic  field  by  requiring  that  the  field 
satisfy  equations  (1.29)  and  (I.30)  on  the  domain  and  that  all 
the  fields  have  the  Identical  two  periods.   The  two  periods 
which  we  specify  are  stated  In  terms  of  flux  conditions;  they 
are  the  flux  0  and  the  value  H  of  the  flux  function  on  the 
fixed  boundary.   Thus,  If  B  denotes  the  magnetic  field  on  the 
unvaried  domain  D,  with  flux  function  ^,  then  we  must  have 

(15)  0,  =JJ[ ^^Jrara^ 

and 

(16)  ^  =  H 

on  the  fixed  boundary  C  of  D.   Similarly,  If  we  represent  the 
magnetic  field  on  the  varied  domain  D  by  B  and  Its  flux 
function  by  ip^ ,    then  the  periods  In  this  case  are 

pn/krTp*  +Bt\ 

(17)  0^  =  //  ^1^  rdrdcf) 


D^ 


a 


and 

(18)  ^*  =  H 


* 


on  the  boundary  C. 

In  addition,  we  also  Introduce  an  energy  for  each  magnetic 
field.   For  the  magnetic  field  B  on  the  domain  D  it  is 
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(19)  ^  ^  ^JJ  B^^'^^^'l^  ' 

D 
while  for  the  magnetic  field  B  defined  on  D   it  is 

(20)  E*  =  -^  //   B*^rdrd(t)  . 

Our  aim  is  to  develop  an  expansion  of  the  energy  E*  in 
terms  of  E  and  in  powers  of  e.  We  will  establish  both  upper 
and  lower  bounds  on  E*  to  show  that  the  first  variation 


(21)  5E  =  e  lim  ^  "  ^ 

e-^0    ^ 


is  equal  to 

(22)  BE  =  r  (^  -  p)6vds 


The  requirement  that  the  appropriate  domain  for  the  solution  of 
the  free  boundary  problem  minimizes  the  magnetic  energy  Implies 
the  condition  that  the  first  variation  5E  vanish  for  all  per- 
turbations 5v.   Therefore  we  must  have 


^  -  P  =  0 


at  each  point  of  the  free  boundary  r. 

First,  we  will  obtain  an  upper  bound  on  E  .   Since  both  B 
and  B*  have  zero  divergence,  it  is  possible  to  introduce  the 
flux  functions  ip   and  f     respectively  for  each.   Expressed  in 
terms  of  these  flux  functions,  the  energies  (19)  and  (20)  are 
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(23) 

and 

(24) 


1 

2  J. 

D 


=*  -  ki 


/2  ,2  ^2' 

li  +  _£  +  _i 

2  ^  ^  ^   2 

r  a  a 


rdrdcf) 


D^ 


\pf        ^*2    B*^' 

2    ^^  ^   2 
r     a     a  i 


rdrd(j) 


Since  the  flux  0  Is  identical  for  both  magnetic  fields,  we  must 
have 


(25) 


5 — ^  rdrdcp 


kr^J  +  B»  '^ 


D 


a 


rdrdif) 


D^ 


We  now  introduce  a  third  function.   It  is  the  comparison 
function  i\i      (r^cf))  defined  on  the  domain  D*  by 


(26) 
Clearly, 


r^{v\k*)    =   ^(r,(t)) 


f**  =   H 


on  the  fixed  boundary  of  D  and 


^**   =   0 


on  the  free  boundary  of  D'*,   In  order  to  have  the  flux  0     of  ^** 
identical  with  that  of  the  other  two  functions,  we  define  the 
constant  B^*  so  that 


C 


(27) 


0   = 


D' 


kr^^  +  B^ 


a 


rdrdcp  . 
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Furthermore,    the   energy  of   this   magnetic   field   Is   defined  as 


(28) 


^** 


+ 


2+2 

o  o 


rdrdii    . 


Since  tp        Is  defined  on  the  same  domain  as  f*   and  has  Identical 
periods,  the  Dlrlchlet  principle  of  Section  2.2  can  be  applied 
to  show  that 


(29) 


E*  <  E** 


*■)(■ 


Next,  we  proceed  to  express  E   In  terms  of  E.   We  start 


with 


E 


*-Sf      1 


D 


■X-   I- 


**2     ^^2  ^^2 

4-  i +  — 2 

*2      ^2      *2 
r*     a  o 


■K  -fr  c~  i 


r*dr*d({)*  , 


where 


o*^  =   1  +  k^r*^ 


This  Integral  can  be  transformed  by  a  change  of  variables  Into 


E 


^* 


2, J 


D 


r,2 

,2     7 

^^ 

+ 

^r* 

*2 

*2 

|_r 

c^      J 

B 


■x-* 


Jr*drd(|)  +  -^  //  -§2  r^dr^dc})*  , 


where 


Mr%r) 


S(r, 


) 


From  equations  (10)  and  (11),  we  find  that 
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5r_ 


hs. 


=  1  -  e 


^+0(s^)  , 


hh 


hs. 


and 


S] 

*:   -    -    -^ + 

5< 

+   O(e^) 


Direct  calculations  then  yield 


/as       as  \ 


and 


i^ 


Sr    Scj) 


as- 


as^  \ 


^^^^-M^rar"^^af;"°^^^ 


^'i    ,    .,     ^^2^,  .,_2. 


V  =  ^r-^[^r^  +  V'^^  +0(e^) 


In  addition,  since  the  flux  is  constant,  algebraic  manipulation 
shows  that 


*-D*^ 


r*B 


dr^dcf)*  = 


/  -^  drd(}) 


D' 


D 


-  ='J7  = 

D 


j^  as^      ^  as^      (1-  k^r^)s- 
"2  ^^  +  -2  s^  + n = 


o 


drdff) 


+  2  7j  B^k 

D 


''^  SS^    r^  aS2    2r^' 


2  ^^  -    2  ^ 


+ 


a 


a 


drd({) 


Hence,  we  obtain 
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(30)   E^ 


D 


r  2     2    2  ' 

a     a 


2^r^<l 


^S, 


^V^    _    ^^    ^C^ct) 


a 


+ 


2      2     2        2 

T/zT    r^    B^r    2kr  ?//  B^ 

r    —?"    2        2 


+  S, 


^f    (1-  k^r^)^^    (1-  k^r^)B^ 


^ 


+ 


a 


a 


C  r 

— n — 


drd(t)  +  0(e) 


Next,  we  apply  Green's  theorem  to  equation  (50)  in  order 
to  express  it  as  a  line  integral.   Remembering  that 


%  +  H^fl  -^^^-0 , 


this  transformation  yields  the  identity 


(31)      E 


** 


/2  ,2         „2 


\ 


Sd  ^ 


dcf) 


+   S- 


/^Vi 


dr  +  S, 


2r^^,         2B^kr  ^ ,  \ 
o^  +  ^ ^     d({) 


a 


2  2  2  2 

f\         v-ip  2B^kr  t//  B^r  .        ^  ^ 


a 


o        ' 


Happily,  this  expression  can  be  simplified.   Since 


^s  =  ° 
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on  the  boundary  of  D,  we  must  have 

Also,  the  fact  that  the  tangent  vector  and  the  Inner  normal  are 
orthogonal  implies  the  relationships 


and 


r  .  =  -ri 
V     ^s 


r  =   ri 
s     ^v 


Using   these   inequalities,    we   can  bring   equation    {j>l)    into   the 
simpler   form 

(32)        E«*   =  E  +  /(s,  ^  .   r%  mi%  .  %  +  !|lds   +  0(e2) 


However,  further  simplification  is  possible.   The  first  factor 
under  the  integral  is  equal  to  6v;  the  second,  to  -p- .  Further- 
more, the  integrand  is  zero  on  the  fixed  boundary  C  of  D  since 


S^  =  S2  =  0 


there.   Therefore,  it  follows  that 

(33)  E*""  =  E  +  J    ^  6vds  +  0(£^)  . 

r 

If  we  now  combine  this  result  with  equation  (29),  we  find  that 
(3^)  e""  <  E  +  J    ^  5vds  +  O(e^)  , 

r 

the  desired  upper  bound  for  E"^. 
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2A     The  Derivation  of  a  Lower  Bound  for  the  First  Variation 

The  same  methods  which  we  employed  in  Section  2.3  can 
also  be  applied  to  derive  a  lower  bound  for  E  .   However,  it 
Is  first  necessary  to  make  several  additional  definitions.   In 
Section  1.3?  we  defined  the  current  of  the  magnetic  field  in 
the  {r,d)    plane  by 

(35)  I^  =/(B^dr+rBgde), 

where  ^  is  an  arbitrary  closed  curve  lying  in  the  domain  D. 
This  definition  may  also  be  expressed  in  terms  of  the  flux  func- 
tion tp ,    where  it  becomes 

n/  ^1       VTp                 kr  Bp   \ 
(56)        I^  =  /f  ^dr  -  ^d^  +  2^dr   . 

Furthermore,  the  value  of  this  integral  is  independent  of  the 
path  chosen  since  ip   satisfies  equation  (1.35)- 

Next  we  define  the  constant  C^  and  the  function  xi^A)    '^^ 
the  domain  D  by 

B. 

(37)  ^r  =  T^ 
and 

(38)  X.  = 


r   rl 


2 

VTp  kr  C^ 

(39)  ^  =  --^^ — ^ 
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A  comparison  of  the  definition  of  x   with  equations  (1.48)  and 
(1.49)  shows  that  x  is  the  potential  function  O,  normalized  so 
that  it  has  current  one.   That  is. 


(to) 


z 


We  also   define   the   energy  M  of   x  ^s 


(ti) 


M  =  i 


r 

D       -  ^ 


rdrdcj) 


Moreover,  making  use  of  equation  (40),  it  follows  that 


(^2) 


M  = 


E 

7^  • 


It  is  an  easy  matter  to  make  similar  definitions  for  the 
magnetic  field  B*  defined  on  the  domain  D*.  In  this  case,  the 
current  is 


(45) 


I*  =  /   B*dr  +  vB*a9    , 
z       J         V  9        ' 


where  H*   is  a  closed  curve  lying  in  D*.   Expressed  in  terms  of 
the  flux  function  ip* ,    this  definition  becomes 


(44) 


I*  = 


^1      r^ 


^  d(l)  +  — ^  d(t) 


H    ■ 


*  . 


Continuing  in  a  parallel  fashion,  the  constant  C^  is  defined  by 


B  y 
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and  the  normalized  potential  function  x     ^Y 


and 


^ 


y*  =  -^ 


o 


*        -! 


In  addition,    the   energy  M      is   defined  as 


(45)  M^   =.  1 

2 


■K-c 

xf  +%+    (C^-kxp2]rdrd(t)    . 


D^ 


If  n     denotes   the   potential   function  for  B*   then,    clearly. 


X*(r,<.)    "52%ii 


and 
(46) 


*        E 


M      = 


■*2 


The  flux  0    ,    which  is  identical  for  the  functions  ^  and 
"ii* ,    may  also  be  expressed  in  terms  of  either  of  the  normalized 
potential  functions  x  ^^^   X* '      Therefore, 


0=1 
'^z     z 


(-kx^  +  C^)rdrd(^ 


or 


0  =  I* 
z     z  J 


(-kxJ  +  C^)rdrd^ 


D' 
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Now  that  these  definitions  have  been  made,  we  proceed 
as  we  did  in  the  last  section.   A  comparison  function  x**   is 
defined  on  the  domain  D*  with  energy 

-)f  -)f  £_ 

X^^  +\-+   (C^-kxJ")2]rdrd^  . 


(47)    M*^  =  -^j 


D* 


On  one  hand  an  analogue  of  Dirichlet's  principle  will  allow  us 
to  compare  M**  and  M*,  while,  on  the  other  hand,  we  can  relate 
M**  to  M  by  elementary  calculus.   Thus,  we  obtain  a  relation- 
ship between  M  and  M"^  which,  when  used  in  conjunction  with  equa- 
tions (^2)  and  (46),  enables  us  to  derive  the  desired  lower 
bound  on  E'*^. 

The  comparison  function  x**   is  defined  by  extending  the 
function  x   onto  the  domain  D*.   On  the  subdomain  Dn  D*,  we 
merely  set 

X**{r,^)    =   x(r,^)  . 

The  extension  onto  D*'-  D,  however,  is  not  as  easy.   On  the 
boundary  of  D  we  have 


^     =  0 
^s 


which,  when  expressed  in  terms  of  the  function  x>    becomes 


2 


(^8)  rxj^  -  ^4^  +  krC.r^  =  0 


¥e  extend  x(r,<}))  onto  D*-  D  by  requiring  that  equation  (48)  be:. 
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satisfied  along  each  outer  normal.   Moreover,  the  constant  ct* 


** 


associated  with  x   is  taken  to  be  identical  with  C^. 

From  this  construction  we  know  two  important  facts  about 
X** '      First,  the  function  is  continuous.   Second,  the  current 
is  one.   We  will  now  use  these  facts  to  prove  an  analogue  of 
Dirichlet's  principle  known  as  Thomson's  theorem,  which  relates 
the  energies  M*  and  M**.   Let  the  function  f(r,(t))  be  defined 
on  D*  by 


Since  both  x**   and  x^  have  period  one,  f  must  be  single-valued. 
We  can  then  express  the  energy  M*"^  in  terms  of  f  and  x*^ 
obtaining 

,2 

rdrdd)  . 


M*^  =  \ 


(f,^X;)'.^^^^MC^-lcx|-kf^)' 


Consequently, 


f? 


(49)      M**    =M-    +|/nf2    +4+   k^f^Wdct) 


f      * 


7^  KxJ  +  4^  +  k%%  +  kc^%  Wd^ 


D 


•X- 


+  \II  ^^r  ^T  -^^C^^  -2C^kf^+2C*kx|-2kCjf^)rdrd<t) 


D* 


However,  since 


^s 
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X*  satisfies  the  boundary  condition 


^xl^ 


r^  s 


2  * 


krCjr  =  0 
C  s 


on  both  boundaries  of  D*.   Prom  this  equation  we  derive  the 
integral  relation 


^f  y* 


CH 


rdrdcj)  =  0 


by  an  application  of  Green's  theorem.   We  use  this  identity  to 
simplify  equation  (49)  and  obtain 


(50)   M 


** 


f2 
4^4^    ((^C-^C^-^'^)  jrdrd^ 


+  Jf  [  -2C^kxJ  +  2C*kxJ  +  2C^cJ  -  2C^^]rdrd<t) 


D^ 


Equation  (50)  can  be  rewritten  in  terms  of  the  flux  0  as 


(51)    M**  =  M*  + 


(C^-  Cp0^ 


+ 


'^-.li 


D* 


fj  +^  +  ((c^  -C*)  -kf^)' 


rdrd<{) 


We  can  neglect  the  last  term  in  equation  (51)  since  it  is 
positive  and  thereby  obtain  Thomson's  theorem,  that 


(52) 


(C.  -  Ct)0 
M**  >  M*  +  ^ ^ — - 
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We  can  also  evaluate  M   in  terms  of  M.   Since  we  defined 


X        =   X 


on  the  Intersection  of  D  and  D  ,  it  is  possible  to  write 


M 


** 


1 
2 


DnD* 


x'+4+  (c^-^x^)' 


rdrd(j) 


-I 


D*-D 


2    Xj 


^*i 


X**-  + 


+  (C^-kx|*)' 


rdrd(j) 


It   then  follows   that 


**2 


(53)    n^*=n-^^ff    [x^^  +^+  (c^-kx|^)' 


rdrdtj) 


D*-D 


1 
2 


,2.  a 


D-D' 


Xr  +  4  +    ^^C  -kX<},)^|rdrd^ 


These  two  integrands  may  be  evaluated  by  writing  them  in  (v,s) 
coordinates.   Under  this  transformation,  equation  (53)  becomes 


*^2 


M^*   =  M   +  i    ff     \x;*^   +  \-  +    (C^  -  kxf) 


1 
2 


where 


D*-D 


D-D' 


rJdvds 


Xr+^+    (C^-kx^)2 


rJdvds    , 


(5^) 


J  = 


r  r 

V  s 


^v        ^s 


h  +  0(e) 


r 
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The  Inner  integration  may  be  performed  approximately.   Both  x 
and  X**  3-^s  continuous  on  their  respective  domains,  and  we  may 
approximate  their  values  along  the  inner  normals  to  order  e  by 
their  values  on  the  boundary  5D.   Hence, 


M*"^  =  M 


1 
2 


r 


^r  +   2  +  (C^-kx^)^j6vds  +  0(£2) 


and  therefore 


(55) 


M"^*  =  M 


21 


B^5vds  +  O(e^) 


By  combining  the  identity  (55)  with  Thomson's  theorem, 
equation  (53) j  we  establish  that 


(56) 


M 


21 


1    /  ^2 


B^6vds  >  m""  +  —^ ^  +  Oie"^)    . 


It  now  remains  to  relate  M  and  M*  to  E  and  E*.   We  can  write 


E  = 


B^0    -, 
2     2 


"Pk        ^ 


'A 


D 


a 


krB^^^ 


^ —  rdrdcj)  , 


where  0  is  defined  by  equation  (15) •   An  integration  by  parts 
then  yields  the  important  identity 


(57) 


B;„0      HI 
E  -  ^—  +  -^  , 


where  H  is  the  value  of  ?//  on  the  fixed  boundary  and  I   is  the 
current.   Moreover,  a  similar  manipulation  of  the  definition  of 
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E*  yields  the  parallel  result 


(58)  E*=^  +  S^ 


Since  X  is  obtained  from  tp   by  dividing  by  I  ,  we  clearly  have 

l2  -     21    -      4E 

and  similarly 

^*  H  +C^0     (H  +  ct0  )^ 

1*2      21*        4e* 
z         z 

We  now  express  relation  (56)  in  terms  of  E  and  E   to 
obtain 


"^^z  r  z     z 


Furthermore,  this  inequality  can  be  inverted,  yielding 

(59)    ^^-^ ^ ^ 1 ^ ;r  +  °^^^^ 

By  expanding  the  right  side  of  equation  (59)>  we  obtain 

4E      ,      i6e^        r^2.  ..  _  ^^z  ,  ^,  2. 


-p  +  — ^-^^^^^ jr  /  B'^Bvds  < +  O(e^) 

(H+C^0^)^    2I^(H  +  C^0^)^  -'^  H-C*0^  +  2C^02 


Continuing  our  algebraic  manipulation,  we  see  that 
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IJ(H+C.0  )^  1  r  2  2 

(60)     i  ? ^^ E  >  ^  /  B^bvds  +  O(e^)  . 

(H-C*0^+2C^0^)        ^1 


This  inequality  may  be  rewritten  as 

,2 


E*(H  +C^0  ) 


^ o  -  E  >  ^  /  B  6vds  +  0(e) 


(H+C^0^)^-0^(C^  -Cp'  p 


Since 

C^  -  Cp  =  0(e)  , 

2  2 

the  term  (C.,  -  C^)   is  of  order  e  and  we  may  neglect  it.   There- 
fore, 

E*  -  E  >  |-  /  B^5vds  +  O(e^)  . 


Moreover,  a  combination  of  this  inequality  with  equation  {j>6) 
yields  the  result 


(6l)  6E  =  -^  /  B^Svds  . 

r 


The  side  condition  that  the  areas  of  D  and  of  D  must  be 
identical  can  be  expressed  in  a  more  transparent  way.   Let  A* 
denote  the  area  of  D*  and  A  the  area  of  D.   Then 


(62)  A*  =  J  I    r*dr*d(l)* 

=  //  r*Jdrd(t)  , 
D 


^9 


where  J  is  defined  in  equation  (5^)-   Expanding  equation  (62),  we 
find  that 


A   =  A  +  e 


D 


Since 


we  must  have 


A^  =  A  , 


drdcj)  +  0(e)  . 


D 


SS      BSp 


^ 


1 


drdcj)  =  0 


The  divergence  theorem  is  applied  to  this  equality  to  obtain 


/  (rS-j_d(t)  -  rSgdr)  =  0  , 


or 


(65) 


6vds  =  0 


Hence, the  condition  that  the  area  of  D  is  identical  to  that  of  D 
is  equivalent  to  equation  {6^).      We  use  the  method  of  Lagrange 
multipliers  to  add  equation  (63)  as  a  constraint  to  our  varia- 
tional problem.   We  are  thus  led  to  the  necessary  condition 


(64) 


(^  -  A)6vds  =  0 


in  order  for  the  free  boundary  r  to  minimize  the  energy.   However, 
this  is  equivalent  to  the  free  boundary  condition 


B 


at  each  point  on  r. 
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Chapter  3  -  Conformal  Mapping 

3. 1   Introduction 

Although  a  large  body  of  knowledge  exists  regarding  the 
solution  of  a  partial  differential  equation  on  a  fixed  domain, 
much  less  is  known  when  a  free  boundary  condition  is  introduced. 
The  basic  difficulty  in  developing  a  numerical  scheme  is  to  be 
able  to  handle  the  changing  shape  and  location  of  the  free 
boundary.   In  the  method  to  be  described,  we  introduce  an 
auxiliary  domain  R  with  fixed  boundaries  in  the  (u,v)  plane  and 
use  the  method  of  steepest  descent  to  find  a  mapping 

X  =  x(u, v) 

y  =  y(u,v) 

which  maps  R  conformally  onto  an  arbitrary  domain  in  the  (x,y) 
plane.   The  solution  to  the  free  boundary  problem  is  then  found 
in  a  parametric  form  which  depends  on  this  conformal  mapping. 
Thus,  the  problem  of  computing  the  location  of  the  free  boundary 
is  transferred  to  that  of  finding  such  a  conformal  mapping. 
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3- 2  The  Dirichlet-Douglas  Integral  and 
the  Method  of  Steepest  Descent 

Let  us  consider  a  free  boundary  problem  which  is  a 
generalization  of  the  one  described  in  the  first  chapter.   As 
there,  we  let  D  denote  a  domain  with  fixed  boundary  C  and  free 
boundary  r.   However,  we  generalize  the  partial  differential 
equation  to  be  solved  on  D  to 

(1)  ^^^^y^^xx  +  2b(x,y)?^^y+  c{x,y)V/yy  +d(x,y)^^ 

+  e(x,y)^y  +  f (x,y)?^  =  g(x,y)  , 

where  the  functions  a,  b,  c,    d,  e  and  f  are  all  sufficiently 
smooth  in  both  x  and  y.   The  differential  operator  is  assumed 
to  be  elliptic,  that  is, 

a(x,y)-c(x,y)  -  b  (x,y)  >   0 

on  the  domain  D.   If  we  denote  the  above  elliptic  operator  by 
L,  then  the  general  free  boundary  problem  is  to  solve 

(2)  L(V/)  =  g(x,y) 

on  the  interior  of  D,  together  with  the  boundary  conditions 

(3)  ^  =  ^1 
on  C  and 

(4)  Tp  ^  0 
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on  r.   An  appropriate  free  boundary  condition  is  also  assumed 
to  hold  on  r,  but  we  shall  not  specify  it  here. 

As  is  well  known,  the  elliptic  operator  L  can  be  trans- 
formed into  the  canonical  form 


(5)  L(^)  =  A^  +D^^  +E^^  +  f^ 


by  the  introduction  of  new  independent  variables  u  and  v.   The 
requirement  that  the  transformation  reduce  L  to  canonical  form 
Implies  that  the  functions  u(x,y)  and  v(x,y)  must  satisfy  the 
pair  of  first  order  linear  equations 

bv  +  cv 

(6)  u  ^    ^     y 

/ac  -b 
and 

av  +  bv 

(7)  u  =  -    ^    _y  . 

j  ac  -  b 

This  pair  of  first  order  equations  for  u  and  v  is  known  as  the 
Beltrami  equations. 

Since  it  will  prove  to  be  more  convenient  to  regard  u  and 
V  as  the  independent  variables,  we  interchange  the  independent 
and  dependent  variables  in  equations  (6)  and  (7)  to  obtain  the 
inverse  Beltrami  equations 


-  ay  +  bx 
(8)  --       ^    ^ 


X 

V 


,2 

ac  -  b 


and 
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ex  -  by 

(9)  y,  =   "   " 


^ac  -b' 


The  domain  R  of  these  two  equations  is  the  image  of  D  under 
equations  (6)  and  (Y)* 

To  simplify  matters,  let  us  introduce  the  complex  variables 

z  =  X  +  iy 

and 

w  =  u  +  iv  . 


Complex  differential  operators  -v—  and  -—  will  be  introduced  by 

^     Sw 
the  relations 


and 


5  _  1  f  S      S  ^ 
dw 


Using  this  notation,  we  can  rewrite  the  pair  of  first  order 
equations  (8)  and  (9)  as  one  second  order  equation. 


2  2 

ay  -  2by  x  +  ex  ,  =  0 
'^w    ''w  w    w 


Moreover,  if  we  define  the  function  f(z)  by 


(10)  f(z)  ==  ay^  -  2by^x^  +  cx2  , 


then  the  problem  of  solving  the  inverse  Beltrami  equations  is 
equivalent  to  finding  the  transformation  z  =  z(w)  which  makes 
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(11)  f(z)  =  0 

on  the  domain  R. 

The  complex-valued  function  z(w)  has  the  additional 
property  that  it  is  a  conformal  mapping  of  the  domain  R  onto 
the  domain  D  under  the  metric 


2  2  2 

ds  =  a(x,y)dy  - 2b(x,y)dxdy  + c(x,y)dx 


We  make  use  of  this  fact  and  cite  the  following  theorem  of 
Morrey  [17].   Let  O  be  the  class  of  dlf ferentiable  functions 
z(w)  which  map  the  domain  R  in  the  w-plane  onto  a  given  domain 
D  In  the  z -plane  and  map  the  boundary  of  R  onto  the  boundary 
of  D.   The  desired  conformal  mapping  from  R  onto  D  is  that  func- 
tion which  minimizes  the  Dlrlchlet-Douglas  Integral 

(12)      ^  =  "I  / /  [a(Vy)^-  2bVy.Vx  +c(Vx)^]dudv  . 

R 

Note  that  we  have  not  specified  the  image  of  the  boundary 
of  R  point  by  point  under  the  mapping  z(w).   Thus,  the  Image  is 
free  to  slide  along  the  boundary  of  D  as  we  move  from  one  func- 
tion in  C  to  another.   It  is  this  degree  of  freedom  which  allows 
us  to  derive  the  Beltrami  equations  as  well  as  the  usual  Euler 
equations . 

Some  topological  considerations  are  pertinent  at  this  time. 
First,  the  order  of  connectivity  of  the  domains  D  and  R  must  be 
identical.   Next,  the  number  of  points  whose  Images  must  be 
specified  to  Insure  the  uniqueness  of  the  mapping  is  also 
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determined  by  the  order  of  connectivity.   For  example,  if  D  is 
simply-connected,  the  images  of  three  points  must  be  specified 
to  insure  uniqueness;  while  if  D  is  doubly-connected,  only  the 
image  of  one  point  must  be  specified.   Last,  if  the  order  of 
connectivity  of  R  and  D  is  greater  than  one,  the  Riemann  moduli 
of  the  two  domains  R  and  D  must  be  identical.   Clearly,  the 
Riemann  moduli  of  R  will  affect  the  location  of  the  free  boundary 
of  D. 

We  now  apply  the  method  of  steepest  descent  in  order  to 
find  the  conformal  mapping  z(w).   To  do  this,  we  introduce  the 
time  t  as  a  parameter  and  suppose  that  x  and  y  depend  on  both  t 
and  w,  so  that 

X  =  x(w, t) 

and 

y  =  y(w,t)  . 

Thus, I  is  a  function  of  the  time  t  that   we  want  to  minimize  as 
t  goes  to  infinity.   This  can  be  done  by  making  I  as  small  as 
possible.   Taking  the  derivative  with  respect  to  time,  we  see 
that 

i   =  ^  {a(Vy)^+  2aVy.Vy  -2bVyVx-  2bVy.Vx 


R 


2b Vy .  Vx  +  c  ( Vx )  ^  +  2c Vx  •  Vx}  dudv 


The  time  derivatives  of  the  coefficients  a,  b  and  c  may  be 
expressed  in  terms  of  the  derivatives  x  and  y  by  the  relations 
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and 


a  =  a^x  +a  y  , 


b  =  b  X  +b  y 
X    Y 


°  "  V  +  V^ 


An  application  of  Green's  theorem  then  yields 

(13)  i   =   -  ff   (y[V(aVy)-  V{bVx)-  |(ay(Vy)^-  2byVy •  Vx  +  c^ ( Vx)^)  ] 
R 

+  i[V(cVx)  -  V(bVy)-  -l{a^(Vy)2-  2b^Vy.Vx+  c^(  Vx)^)  ]  }dudv 


(y(ay^  -bx^) +x(cx^  -  by^)}ds  , 


bR 


where  s  denotes  the  arc  length  along  SR  and  v  denotes  the  Inner 
normal. 

The  paths  of  steepest  descent  are  those  which  make  the 
integrands  on  the  right  positive.   On  the  interior  of  R,  this  is 
accomplished  by  setting 

X  =  V(cVx)  -  V(bVy)  -  1  [a^(Vy)^  -  2b^Vy.Vx+  c^(Vx)^] 


and 


y  =  V(aVy)  -  V(bVx)  -  |  [ay(Vy)^  -  2byVy.Vx  +Cy(Vx)^] 


Since  the  method  of  steepest  descent  tells  only  the  "direction" 
to  move,  but  not  the  amount,  these  formulas  may  be  multiplied  by 
arbitrary  positive  constants  a.,  and  a^   to  obtain 
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(14)  k   =   a^{V(cVx)  -  V(bVy)  -  i  [a^( Vy )2  _ 2b^Vy • Vx +  c^( Vx)^] } 
and 

(15)  y  =  a2{V(aVy)  -  V(bVx)  -  |  [a^i^yf-   2byVy .  Vx  +  c^ ( Vx)^] } 

on  the  interior  of  R. 

It  is  a  more  complicated  procedure  to  derive  the  paths  of 
steepest  descent  on  the  boundary  of  R.   Ideally,  we  would  choose 


and 


x  =  ex  -  by 
V    "^  V 


y  =  ay^  -  bx^ 


on  the  boundary  of  R.   However,  any  iterative  scheme  must  satisfy 
the  requirement  that  the  boundary  of  R  maps  into  the  boundary 
of  D.   Thus, a  point  which  lies  on  the  boundary  of  SD  must  remain 
there,  and  its  motion  must  be  restricted  to  a  path  along  the 
boundary.   This  is  accomplished  by  projecting  the  path  prescribed 
by  the  method  of  steepest  descent  onto  the  vector  (x  ,y  ),  the 
tangent  of  5R.   Thus, the  actual  iterative  scheme  used  is 

,  .     7[ay  y  -bx  y  -  bx  y  +  ex  X  ] 

(16)        (x,y)  =  ^-^ "^ ^ ^^  (x^,y^) 

ay^  +  cx^ 

where  y   is  an  arbitrary  positive  constant. 
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3. 3  The  Steady-State  Solution 

It  Is  still  necessary  to  check  that  the  solution  of  equa- 
tions (l4),  (15)  and  (16)  Is  also  a  solution  of  the  original 
problem.   As  t  -*  00  ,  a  steady  state  is  reached  and  we  have 

(17)  0  =  V(cVx)  -  V(bVy)  -  i  [a^(Vy)^-  2b^Vx.Vy  +  c^(Vx)^] 
and 

(18)  0  =  V[aVy)  -V(bVx)  -  1  [ay(Vy)2  -  2byVy •  Vx  +  c^ ( Vx)^] 

on  the  interior  of  R  and 


(19)  0  =  ay  y  -  bx  y  -  bx  y  +CX  X 

^  -'''  '^v  s     V  s     s'^v     V  s 


on  the  boundary. 

For  convenience,  we  specify  the  domain  R  so  that  its 
boundaries  are  either  horizontal  or  vertical  lines  in  the  (u,v) 
plane.   In  particular,  if  R  has  connectivity  one,  then  we  take 
R  to  be  be  the  rectangle  lying  in  the  first  quadrant  of  the 
(u,v)  plane  with  one  corner  at  the  origin.   If  R  is  of  connec- 
tivity tvjo,  then  we  modify  this  definition  so  that  R  is  a 
finite  strip  lying  on  the  u  axis  and  periodic  in  u.   For  higher 
orders  of  connectivity,  a  more  complicated  figure  must  be  chosen 
for  R. 

Combining  equations  (l?)  and  (18)  yiedls  the  result  that 


(20)  f _  =  0 

z 
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and  therefore  that  f{z)  is  analytic  on  the  domain  R.   Equation 
(19)  is  equivalent  to 

Im  (f)  =  0 

on  the  boundary  of  R.   Since  f  is  analytic,  Im  (f)  must  be  a 
harmonic  function,  and  an  application  of  the  maximum  principle 
shows  that 

(21)  Im  (f)  =  0 

throughout  the  domain  R.   Combining  this  result  with  our 
knowledge  that  f  is  analytic,  we  see  that 

(22)  Re  (f)  =  H  , 

a  constant,  on  the  domain  R.   If  we  could  show  that  H  =  0,  then 
we  would  know  that 

f(z)  =  0 

on  the  domain  R  and  that  z(w)  is  the  desired  conformal  mapping. 
To  do  this  requires  the  use  of  the  free  boundary  condition, 
which  we  will  discuss  in  the  next  chapter. 

It  should  be  noted  that  the  generality  of  the  description 
of  the  coefficients  of  the  operator  L  enables  us  to  make  a 
possible  simplification.   If  h(x,y)  is  an  arbitrary  function, 
continuous,  and  never  zero  on  D,  then  the  free  boundary  problem 
with  the  differential  equation 
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hL(V/)  =  hg 

has  exactly  the  same  solution  as  our  original  problem.   Multi- 
plication by  such  a  function  h(x,y)  may  be  of  advantage  in  some 
way,  such  as  increased  stability  or  quicker  convergence  of  the 
scheme. 
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Chapter  4  -  The  Complete  Algorithm 

k . 1      Introduction 

We  now  proceed  to  specialize  the  algorithm  for  the 
minimization  of  the  Dirichlet-Douglas  integral  (3.12)  to  the 
free  boundary  problem  defined  in  Chapters  1  and  2.   Since  this 
problem  has  already  been  stated  in  terms  of  the  variables  r 
and  (f),  we  will  continue  to  use  them,  and  not  x  and  y.   The 
coefficients  of  the  differential  equation  (3-1)  are  then 


a  =  — 2  ,      b  =  0  , 
a 

1  ,  (1-  k^r^) 

c  =  -g  .       d  =  ^ j^ '- 

V  ro 

e   =  0  ,               f  =  0 


and 


2kB^ 
o 


For  this  set  of  coefficients,  the  inverse  Beltrami  equations 
are 


(1)  ov^  =   r^^ 


(2)  ar^  =  -rcf)^ 


and  the  complex-valued  function  f  (3- 10)  becomes 
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(3) 

where 

and 


f(z)  = 


(b    r 
^w  ,   w 
-J  +  -^  , 

a    r 


z  =  r 


+  icj) 


w  =  u  +  Iv 


The  pair  of  equations  (l)  and  (2)  is  the  analogue  of  the  Cauchy- 
Riemann  equations  under  the  metric 


,2    1   ,i2  ^1,2 
ds  =  —^  d^      +  — ^  dr 

o  r 


It  follows  that  the  iterative  scheme  we  devised  in  Chapter 
3  for  the  minimization  of  the  Dirichlet-Douglas  integral  (3-12) 
simplifies  to: 


and 
(5) 


=  a^(V(J^  Vr)  -  I  [(\)  (Vct))^  -  {\)    (Vr)^]} 


I   =  ctJV{-^  Vc}))} 


on  the  interior  of  R  and 


(6) 


{r,(t)) 


7 


^  K.^c  +  a  r  r 


V'  s 


V  s 


2:2  ,   2~2" 

r  (b  +  a  r 
^s     s 


(-3'*.) 


on  the  boundary  of  R.   Furthermore,  the  fact  that  the  domain  D 
is  doubly-connected  implies  that  the  auxiliary  domain  R  must  also 
be  doubly-connected. 
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When  a  steady  state  solution  of  equations  {h) ,    (5)  and 
(6)  Is  achieved,  we  know  that  the  function  f  has  the  properties 

Im  (f)  =  0 
and 

Re  (f)  =  H  , 

where  H  Is  a  constant.   Unfortunately,  these  two  conditions  are 
not  sufficient  to  show  that  the  mapping  r  +  l(j)  Is  conformal.   It 
Is  also  necessary  that  the  constant  H  be  zero.   The  free 
boundary  condition,  which  we  neglected  In  Chapter  J),    will  be 
used  to  attain  this  goal. 


4. 2  Application  of  the  Method  of  Steepest  Descent 
to  the  Free  Boundary  Variational  Problem 

In  Chapter  2,  we  showed  that  the  free  boundary  condition 
for  the  equilibrium  configuration  of  a  helically  symmetric 
magnetic  field  could  be  expressed  as  the  solution  of  a  varia- 
tional- problem.   This  problem  consisted  in  minimizing  the 
functional 

(7)  "^   "  "I  / /    Isl^rdrdcj)    -   A^    //   rdrdcj) 

D  D 

over  all  domains  D,  where  B  was  a  helically  symmetric  magnetic 


field  satisfying  V-B  =  0  on  D,  VxB  =  0  on  D,  B-v  =  0  on  the 

2 

boundary  of  D  and  A  was  a  Lagrange  multiplier.   In  Sections 

2.3  and  2.4  we  expressed  B  in  terms  of  the  flux  function  f   and 
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showed  that  the  first  variation  was 


(8) 


5J 


P  L  ^  ^   f7   o 


-   A' 


Bvds 


where  r  is  the  free  boundary  and  5v  is  the  variation  of  r  along 
its  inner  normal. 

Since  we  again  wish  to  minimize  a  functional,  we  again 
use  the  method  of  steepest  descent.   The  functional  J  is  con- 
sidered to  depend  on  the  artificial  time  t  and  we  attempt  to 
minimize  J  as  t  — »■  oo  .   Taking  the  derivative  with  respect  to 
time,  we  obtain 

,2\ 


J 


r  ^ 


2    2 


b: 


-  ^  -  +  -IL  +  ZL 

2  V  2     2     2 

TOO 


A 


v^ds  , 


where  v,  is  the  rate  of  change  of  the  free  boundary  r  in  the 
direction  of  its  inner  normal.   It  then  follows  that  the  path 
of  steepest  descent  is 


(9) 


ip 


€ 


B 


r    a    a 


where  5  Is  a  positive  weight.   The  convergence  of  equation  (9) 
implies  that 


\r4 


;2    ,2    ^2 
tpiip  B 


r   .  "C  \   .2 


r    a    a 


X 


at  each  point  of  the  free  boundary  r.   Since  the  above  identity 
is  equivalent  to 
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B 


P  , 


we  have  satisfied  the  free  boundary  condition. 


h .  J>     Derivation  of  the  Complete  Algorithm 

The  two  variational  principles,  one  for  the  conformal 
mapping  and  the  other  to  minimize  the  magnetic  energy,  are 
combined  now  to  obtain  the  complete  algorithm.   On  the  interior 
of  R  we  solve  equations  (4)  and  (5)>  which  are  specified  by  the 
minimization  of  the  Dirichlet-Douglas  integral,  as  well  as 


(10) 


2    2 


lh         +71/ 

^uu      '^vv 


k  r 


3      ^  ^u^v      ^v^u'^ 


2krB^J 

3~ 


the  partial  differential  equation  for  ip   in  canonical  form.   Here 
J  is  the  Jacobian  for  the  change  of  variables. 


J  =  r  i  -  r  h 
u^v   v^u 


The  boundary  values  for  ^  on  R  are  identical  with  the  ones 
on  D.   However,  the  boundary  conditions  for  r  and  \)   are  more 
complicated.   On  the  upper  boundary  of  R  we  solve 


(11) 


irA)    =  P- 


r  (1)  (ti  +  a  r  r 
^v^s     V  s 


2i  2  ^  2  2 

r  (|)^  +  a  r^ 


(^s'^s)  ' 


while  on  the  lower  boundary  of  R  we  combine  this  with  equation 
(9)  derived  from  the  energy  minimization  problem  to  obtain 
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(12)     {r,<t))  =  P, 


r  2,  I   ,  2 
r  (b  (b  +  a  r  r 
^V^s     V  s 


2i  2  ^22 
r  (b  +  a  r 
^    ^s      s 


(-s'^s) 


-  P) 


^ 


V 


B, 


2  2  ,  2,2 
a  r^  +r  .J)^ 


+ 


A 


(%^^v) 


Since  the  function  ip   is  now  assumed  to  depend  on  u  and  v,  we 

2 

have  written  B   in  terms  of  these  variables 


2         2 

^   "   2  2  ,   272  ^   2 

a  r„  +  r  i  ^    a 


This  quantity  has  been  simplified  by  using  relations  (l)  and  (2) 
and  the  fact  that 


f. 


0 


on  the  boundary.   Alternatively,  we  could  have  written 


(13) 


B   = 


_v____ 

2  2  ,   272 

o  r  +  r  (b 
u    ^u 


B 


+ 


c 

2  ' 


using  relations  (l)  and  (2). 

Two  modifications  to  these  equations  remain.   The  first  is 
necessary  to  fulfill  the  claim  that  the  conformal  mapping  z(w) 
sends  the  boundary  of  R  into  the  boundary  of  D.   While  we  use 
this  requirement  in  the  derivation  of  equation  (11),  nothing 
so  far  in  the  numerical  scheme  we  have  proposed  ensures  that  this 
condition  will  be  satisfied.   We  therefore  modify  equation  (11) 
to 
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(U)   (r,^)  =   p. 


r    2i 

r   o   ()     + 

2          1 

2l2  ^ 
r   ()     + 

2    2 
0  r^     J 

(=^s'*s' 


F  +  FT       ^ 
r    9 


where  the  fixed  boundary  of  D  is  defined  by 


(15) 


F(r,(i))  =  0 


The  new  term  in  equation  (l4)  is  another  application  of  the 
method  of  steepest  descent;  this  time  to  ensure  that  the  image 
of  each  point  on  the  upper  boundary  of  R  lies  on  the  fixed 
boundary  of  D.   The  coefficient  pp  which  appears  In  equation  (l4) 
is  a  weight  function. 

Second,  for  numerical  reasons,  we  change  equation  (12)  to 
the  form 


(16)     (r,<t))  =  P 


3 


r  2,    1     , 
r   ()   o     + 
^v^s 

2           1 
^   Vs 

2i2  , 

r  ()    + 

2    2 

a  r^      J 

(-s'^s) 


+  P) 


A^a'-(r  (t)J  +a^rp 


2 ,2  ,  ,  2i2  ^  2  2>„ 


V 


V 


v 


(^.^) 


4.4   Properties  of  the  Steady-State  Solution 

We  still  must  discuss  the  results  when  a  steady-state 
solution  is  obtained  and  show  that  all  the  claims  we  have  made 
are  satisfied.   The  convergence  of  equations  (4),  (5),  (10),  (l4) 
and  (l6)  implies  that  all  the  residuals  are  zero.   Therefore, 
we  have 
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(17)  V{^  Vr)  -  i  (4)  Mf    -  \   (^)  (Vr)2  =  0  , 

(18)  V(J^  Vet))  =  0  , 

a 

2  2  2kB^J 

(19)  A^  -  "^   {^^k^  -  ^4>^)  -   \- 

a  a 


on  the  interior  of  R, 


(20)  ^\<i^v'^^%^v  =  0  ' 


(21)  F(r,(t))  =  0 


on  the  upper  boundary  of  R  and 


(22)  r^(b  (j)  +  a^r  r  =  0  , 

2        2 

(23)  g  2    g,g  +  -|  -  A^  =  0 

o  r^  +  r  (t)^    a 


on  the  lower  boundary  of  R.   It  has  already  been  shown  that  the 
complex  function  f  defined  by  equation  (3)  is  analytic  and 
satisfies 

Im  (f)  =  0 
and 

Re  (f)  =  H 

on  the  domain  D. 

We  now  proceed  to  show  that 

H  =  0  . 
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Using  the  definition  of  f,  we  see  that 

1  ,i2   i2,  ^  1   ,  2    2^ 

H  =  — o   (<P    -  CP   )   +  — 7T   (r   -    ^      J       ' 

a  r 

The  parameter  A  which  occurs  in  equation  (23)  has  not  yet  been 
specified,  so  we  are  at  liberty  to  define  it  as 


(24) 


1 
L 


t 


V 


B 


2- 


2  2  ,  2i  2 
a  r, ,  +  r  ^ 


+ 


u     '^u 


ds 


where  L  is  the  arc  length  of  the  free  boundary  r.   Insertion  of 
this  definition  of  A   into  equation  (23)  shows  that 


(25) 


,2       ^2 

a  r  +  r  6  ,    a 


2^ 

L 


V 


V 


r 


f. 


V 

2  2  ,  27? 
a  r  +  r  d) 
u    ^u 


+ 


B  CO 
a  ^ 


ds 


at  each  point  of  r.   We  now  cite  the  following  theorem  of 
elementary  calculus: 

If  f(s)  and  g(s)  are  continuous  functions  on  the  line  i 
and  if 

f(s)  y  dt  =  j  g(t)dt  , 


then  there  exists  a  point  s  on   H    such  that 

^      o 


f(s^)  =  g(s^)  . 


This  theorem  may  be  applied  to  equation  (25)  to  show  that  there 
exists  a  point  (u  ,v    )    on  the  boundary  r  where 
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;2        ^2        ,2       ^2 

rpfi")       _^    +  -^  -   _  V    +  C 

^^o;         2  2^  2:2  ^   2  ~   2  2  ,   272   ~? 


Hence,  it  follows  that  at  this  point 


fovA  2  2^  2i  2     22^  2,2 

(27)  a  r^  +  r  (j)^  =  a  r^  +r  (j)^  , 


and  consequently  the  value  of  the  constant  H  is  zero.   Thus, 
when  a  steady-state  solution  is  achieved: 


1.   The  function 


^2    2 

f(z)  =  — 2"  "^  ~^ 
a    r 


is  Identically  zero.   Therefore,  the  mapping  z(w)  is  conformal 
under  the  metric 

,2    1   ,.2  ,1,2 
ds   =  — ^  d(j)   +  — ^  dr   , 

o  r 

and  the  inverse  Beltrami  equations 

2      i 

and 

a  r  =  -ri 
V      ^u 

are  satisfied. 

2.   The  mapping  z(w)  sends  the  boundaries  of  R  onto  the 
boundaries  of  D. 
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3-   The  equation  for  ^  in  terms  of  the  independent  variables 
u  and  V  has  the  canonical  form 

1^2   2  2krB.J 

■>p         +TJJ           -    ii^     (^     (J)      -   ^     ({)  )     =    J^ 

a  a 

and  this  equation  has  been  solved  for  tp. 
4.   The  free  boundary  condition 


^  =  p 


is  satisfied  at  every  point  of  the  free  boundary  r. 

Therefore  the  solution  to  the  complete  free  boundary  problem  has 

been  obtained. 


4. 5  A  Simplification  of  the  Algorithm 

The  numerical  solution  of  equations  (4)  and  (5)  is  still 
quite  difficult,  even  though  they  are  to  be  solved  on  a  fixed 
domain.   Both  are  non-linear  and  they  are  coupled  by  the  first 
order  terms  of  equation  (4).   Clearly  it  would  be  advantageous 
to  lessen  these  difficulties,  and  we  do  this  by  means  of  the 
method  described  at  the  end  of  Chapter  J>.      Instead  of  solving 
the  partial  differential  equation 

^r   n^    (l-k^r^)^^   2kB^ 

-2-  +  2  +  4 =  —^     ' 

a  V  va  a 

we  consider  the  equation 
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(28)     *^^  +  (^  +1^^)*^^  +  ^^  =  -^ 

r       ^^       ra         a 


This  changes  the  complex  function  f  to 


2 
(29)  f  (z)  =  ([)2  +  ^  r^ 

r 


and  the  paths  of  steepest  descent  on  the  Interior  of  R  become 


0^   ^\         lfa^\.^.2. 


(30)  r  =  a^[V   °2  ^rj  -  i[  ^  U^^r^ 

\  p    /      ^  r  .r 

and 

(31)  |)  =  a^iy[s/^)]    . 

On  the  other  hand,  both  the  fixed  boundary  conditions  (l4)  and 
the  free  boundary  conditions  (l6)  remain  the  same. 

Comparing  equations  (50)  and  (31)  with  equations  (4)  and 
(5),    it  Is  clear  that  several  advantages  have  been  realized. 
First,  and  most  Important,  the  two  equations  are  no  longer 
coupled.   Second,  the  equations  are  less  complicated.   We  even 
have  the  happy  case  where  equation  (51)  Is  linear. 

When  a  steady-state  solution  of  equations  (28),  (30),  (31), 
(l4)  and  (l6)  Is  attained,  we  repeat  the  proof  of  Section  4.4 
and  observe  the  following: 

1.   A  linear  combination  of  equations  (30)  and  (31)  shows  that 


f   =  0  . 


Hence  f  Is  analytic  on  the  domain  D, 
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2.   Equation   (l4)  implies  that 

Im  (?)  =  0 

on  the  boundary  of  D.   Since  Im  (f)  Is  a  harmonic  function,  it 
follows  that 

Im  (f)  =  0 

on  D  by  an  application  of  the  maximum  principle. 

3-   The  real  part  of  f,  which  Is  the  harmonic  conjugate  of  the 
Imaginary  part,  must  have  the  form 


Re  (f)  =  H  , 


<\ 


where  H  is  a  constant. 

4.   The  convergence  of  equation  {l6)  yields  the  fact  that 


H  =  0 


Therefore  the  function  r+  lif)  is  a  conformal  mapping  of  the 
domain  R  onto  the  domain  D,  and  we  can  use  the  simpler  set  of 
equations  (28),  (30),  (31),  (l^)  and  (l6). 


4.6   The  Finite  Difference  Scheme 

In  order  to  solve  equations  (10),  (30),  (31),  (l4)  and  (l6), 
numerically  we  use  finite  difference  methods  to  approximate  the 
first  and  second  derivatives  of  the  functions  ip,    r,    and  cj).   The 
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domain  R  is  taken  to  be  a  periodic  strip  of  length  A  and 
width  B.   We  impose  on  R  a  lattice  with  uniform  mesh  width  h. 
The  continuous  function  (j)(u,v,t)  is  then  approximated  on  this 
mesh  by  the  discrete  function  h\   ,  ,  where 


(32) 


'^J^k  =  MJh,kh,nAt) 


and  At  is  the  time  interval.   If  (jh,kh)  is  a  point  on  the 
interior  of  R,  we  approximate  the  first  derivatives  of  ^   at  this 
point  by 


^ 


u 


V 


_  VllIliVi^ 


(jh,kh) 


(jh,kh) 


2h 


^j^k+l-'t^j^k-l 


2h 


+  O(h^) 


+  O(h^) 


and  the  second  derivatives  by 


uu 


vv 


(jh,kh)  h"^ 


_  ^J,k+l-^^j,k+^j,k-l 


(jh,kh) 


+  O(h^) 


h 


In  addition  we  approximate  ^   by 


(Jh,kh) 


=   J>^  ^ iiii  +  0(h) 


These  finite  difference  approximations  are  now  inserted 
into  equation  (31)*   Moreover,  we  set 
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At 


and  obtain 


(3.)   ^5.-^) 


(i)l""/_  +  a. 


(K^ri)    +  ^(n)    +  i(n)      ,  (n) 
^  - 


A  further  Improvement  is  possible  on  equation  {33).      We  u 


se 


(34)    ^(-^)  =  ^(^)  .  a. 


r^(n+l)  ^  ,(n) 


^j-l,k  ""   'P,1,k+1  ^  'P.1,k-1 


which  is  the  method  of  successive  point  over- relaxation  for  the 

solution  of  Laplace's  equation  with  relaxation  factor  Qp  >   1. 

On  the  other  hand,  if  a^  <  1  equation  (34)  becomes  a  successive 

point  under- relaxation  method. 

,(n) 


(35) 


A  similar  discrete  function  r\   I   can  be  defined  by 

j^i   =  r(ju,kv,nAt) 


and  the  first  and  second  derivatives  can  be  defined  as  they  were 
for  ^.      The  finite-difference  analogue  of  equation  (30)  is  then 


(36)   r^.'^+l) 


where 


(l-a^)r(-).a^^. 


,^(n+l)     (n+1)     (n)       (n) 


.c    (n)  (n) 
lor;  / o .  / 


(  (n+1)    (n)   ,2   ,  (n+1)    (n)  y 
^"^j+l^k   ''^j-l,k^  +  ^^j,k+l  -^j,k-l'' 


j.k"j,k 


a^.-f  ^l  +  k^r^-f  • 
J^k  j,k 
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To  solve  equation  (10)  for  f ,    we  use  standard  over  relaxation 
methods  and  obtain 


(37)   ^ 


(i^+l)  =  (i_a  l^^'^^  +a  - 

2 


(V/ 


(n+1)  ^  ,,,(n+l)  ^  ^,(n) 


(n) 


j+i,k  +  ^^k+l  -^  ^j-i,k  +  ^-.k-i 


l6a(-fL 


^^j+l,k   ^j-l,k^^^  J,k+1   *Pj,k-l'' 


^^j,k+l   ^j,k-l''^'Pj+l,k   "Pj-l^k^ 


where 


J,k   T(n) 
r    -iS   j,k 


j(n)  _  (^(n+1)  _   (n)   wi(n+l)  _  i{n)       ^ 
'Jj,k  -  ^""j+l^k   ^j-l,k''^<''j,k+l    'Pj,k-1'' 

-  (r^'^+^^  -  r^'^)   )(<i)^'^+^^  -  h^""^      ) 
^^j,k+l   ^j,k-l^^'Pj+l,k   '^j-l,k-'  • 

Both  equations  (1^)  and  (l6)  require  that  we  compute 

normal  and  line  derivatives  on  the  boundary  of  R.  The  following 

finite-difference  approximations  are  used.   If  the  point  (jh,kh) 
lies  on  the  lower  boundary  of  R  then 


(38) 


and 


(39) 


(Jh,kh) 


'j+l,k  -  ^j-l,k  ^  Q(^2j 


(jh,kh) 


2h 


-^^i,k^Hj,,,i-ij,,,,  ^  ^^^2^  ^ 


while  if  (jhjkh)  lies  on  the  upper  boundary  of  R 
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(40)    II        .  ii^l.iS_i±l.ii  +  0(h2) 

^^  (jh,kh)  2h 


and 

=     ^^ ^^^ i^^^  +  0{h^)   , 

(jh,kh)  2^ 


(41)    |i 


Similar  formulas  hold  for  the  line  and  normal  derivatives  of  r 
and  i/   on  the  boundaries. 

On  the  boundaries  we  use  these  finite-difference 
approximations  and  set 

At  =  h  . 
Equation  (l4)  then  becomes 

F    +Ff  ^ 

r   9 

where  T"^  denotes  the  finite-difference  implementation  of 
c 

a  ^J>^  +  o   r  r  ^^ 

p— P p— p —  at  the  j    point  on  the  fixed  boundary  and  where 

r   h     +  a  r 

'  o  c 


s     s 

r  and  <})   are  approximated  by  the  appropriate  form  of  equation 

s       s 

(40).   The  function  F  and  its  partial  derivatives  F   and  Fi  are 
assumed  to  be  specified  independently.   Equation  (15)  is  handled 
similarly,  becoming 
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(t5) 


f-i?^"'^!-:^"'  ^  ^rl^^^i?^-pri^(-3'*.) 


+  Pi^h(A^S^  -  l){r^,^^)    , 


where  T^  denotes  the  finite-difference  implementation  of 

o   't>v^s  +  o   r  r^         ^^ 

^— ^ p— ^ —  at  the  j   point  on  the  free  boundary  and  where 

r  i  +  a  r^ 
^s     s 

r  and  (()  are  approximated  by  equation  (38),  r  and  <{)  by 

So  V         V 

equation  (39) •      Moreover,  Sp  represents  the  finite-difference 
approximation  of 

2,  2i 2  ^2  2, 

o  (r  -j)^  +  a  r^) 

a  ^^  +  (r  (t)^  +0  r^)B^ 


using  equations  (58)  and  (39)  and  A   is  evaluated  by  the  finite 
difference  scheme 


A   = 


N-1  . 


J--L  a 


c^ 


^^^u  +  ^u 


2  +  B^ 


where  N  is  the  number  of  mesh  points  along  the  u-axis  and  i/     is 
evaluated  using  equation  (39) >  r  and  ^     using  equation  (38). 


^.7   Specification  of  Periods 

The  two  periods  which  are  specified  for  the  numerical 
scheme  are  the  current  B^  and  the  flux  H.  However,  when  we 
derived  the  first  variation  of  the  energy,  we  specified  a 
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different  set  of  periods,  the  fluxes  H  and  0  .   Since  the  proper 
formulation  of  the  variational  problem  depends  on  the  periods 
which  are  given,  It  would  appear  that  this  particular  choice  of 
periods  for  our  numerical  scheme  would  Invalidate  our  results, 
which  Include  the  free  boundary  shift  (9)-   We  can  remove  this 
difficulty  by  making  the  assumption  that  at  each  artificial  time 
step  a  new  flux  <fr         Is  chosen  so  that  the  value  B^  remains 
unchanged.   Thus,  we  have  a  sequence  of  fluxes  {0^  '\    which 
converge  to  a  limit  0  as  the  time  goes  to  Infinity. 

An  alternate  approach,  which  was  not  employed,  would  be  to 
hold  0  constant  throughout  the  calculations  and  compute  a  new 
value  of  B^  after  every  Iteration  by  means  of  equation  (1.37) 
or  (1,38).   The  numerical  problem  would  then  be  In  exactly  the 
same  form  as  the  variational  problem  of  sections  2.3  and  2.4. 


4.8  Normalization  of  the  Solutlo 


n 


We  remarked  In  Chapter  2  that  In  order  to  determine  a 
conformal  mapping  uniquely,  the  Image  of  a  certain  number  of 
points  must  be  specified.   In  the  case  under  consideration,  both 
the  domains  D  and  R  are  doubly-connected  and  we  therefore  must 
specify  the  Image  of  one  point  In  order  to  ensure  the  uniqueness 
of  the  solution.   The  question  arises  as  to  how  this  requirement 
will  be  handled.   There  are  two  possible  approaches.   One  is  to 
ignore  the  uniqueness  condition  completely;  the  other  is  to 
specify  the  value  of  one  function  at  a  particular  point,  say 
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Muq,v^)  =  0  . 


This  latter  approach  is  the  one  we  choose  to  follow.   We  specify 
that  at  the  mesh  point  (j,k)  corresponding  to  the  point  (u  ,y    ) 
we  must  have 

for  all  n. 

To  implement  this,  we  use  the  fact  that  the  solution  is 
translation  invariant.   Let  -(^(Ujv),  r(u,v)  and  (})(u,v)  be  solu- 
tions of  equations  (4.10),  [h A)    and  (4.5)  respectively.   Then 
it  is  also  true  that  •i/'(u+e,v),  r(u+s,v)  and  (f)(u+£,v)  are  solu- 
tions of  these  equations  for  an  arbitrary  constant  e.   In 
practice,  we  compute  e  so  that 


(u  +e,v  )  =  0 


and  then  translate  the  functions  ^,  r  and  (|)  by  e  along  the  u-axis, 
Thus  the  solution  is  normalized  in  the  desired  fashion. 


h . 9  The  Choice  of  Weight  Functions 

One  problem  in  obtaining  a  convergent  scheme  is  the  proper 
choice  of  the  weights  a-,  ,  «„  and  a^  which  occur  in  equations 
(34),  (56)  and  (37)  respectively.   We  would  like  to  choose  them 
so  that  the  algorithm  converges  as  quickly  as  possible,  but  not 
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so  large  that  we  overshoot  the  minimum  point  in  the  multi- 
dimensional surface  described  by  the  discrete  formulations  of 
the  Dirichlet-Douglas  integral  {3«12)  and  the  energy  functional 

(7). 

The  three  equations  {j>h) ,    {j>6)   and  (57)  are  non-linear  or 
have  non-linear  boundary  values  which  vary  in  time.   Very  little 
is  known  about  the  proper  selection  of  over-relaxation  factors 
in  this  case.   However,  let  us  consider  the  simpler  case  of 
solving  Laplace's  equation  on  a  rectangle  with  Dirichlet  boundary 
conditions  for  the  discrete  function  x.  ,  .   Applying  the  method 
of  successive  over-relaxation  to  this  problem  yields  the  itera- 
tive formula 


(n+1)         /,       ^     (n)    , 


^.^.I'i^K  A^tl]  +x(^)    ^^x^"^)    n  -^-^."2 
j+l,k        j,k+l        j-l,k        j,k-l  j,k 


(n) 
where  x\  /.  denotes  the  approximate  solution  obtained  after  n 

iterations  and  co  the  relaxation  factor. 

Fourier  analysis  [11]  can  be  used  to  show  that  a  good 

choice  for  o)  is  given  by 


CD  = 


1+  Ch 


where  h  is  the  mesh  width  and  C  is  an  appropriately  chosen 
constant.   With  this  choice  of  cd,  the  rate  of  convergence  of  the 
approximate  solution  x.  /  to  the  steady-state  solution  x\   ,  ^  is 
given  by 


(^^) 


(OD  )       (n)  I   ^  n   - 

X  \  ,   -  X  \  /   <  A  e 
J,k     3,^     - 


Anh 
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where  A  and  A  are  positive  constants  depending  only  on  C. 

Therefore,  the  number  of  Iterations  required  to  achieve  a 

specific  accuracy  grows  like  t-, provided  that  the  relaxation 

factor  o)  is  computed  according  to  equation  (45).   Moreover,  the 

computation  time  T  grows  like  — ^  and  the  total  computing  time 

h^ 
is  given  approximately  by 

(45)  T  =  KMN  , 

where  K  is  the  time  required  to  compute  at  one  mesh  point,  M  is 
the  number  of  interior  mesh  points  and  N  is  the  number  of  itera- 
tions required  for  convergence. 

While  the  theory  cited  is  only  known  to  be  true  in  the 
linear  case,  we  will  apply  it  to  out  non-linear  equations  and  set 

(46)  a 


i    1  +  C.h  ' 

for  i  =  1,2,3.   The  C.  are  appropriately  chosen  constants.   The 
numerical  results  obtained  will  be  used  to  determine  whether 
relations  (44)  and  (45)  are  valid  in  this  case,  and  whether  the 
definition  of  the  a.  given  by  equation  (46)  is  reasonable. 


4.10  The  Initial  Mapping 

The  initial  values  of  the  functions  r.  ,  ,  i  .  ,  and  ij .   , 
are  chosen  as  close  to  the  final  solution  as  possible,  without 
undue  complication.   We  assume  that  the  outer  boundary  of  D  is 
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a  circle  and  assign  the  values  of  r .  ,  and  (j) .  ,  on  the  outer 

J  >  "^       3  }  '^ 

boundary  of  R  accordingly.   Therefore,  the  value  of  r .  ,  is  a 
constant  and  ^  .   -.     is  a  linear  function  in  u.   Recalling  the 
inverse  Beltrami  equations 

o\  =    H^ 


which  the  solution  must  satisfy,  we  integrate  this  system  along 
the  line 

u  =  constant 

to  obtain  the  remainder  of  the  initial  values  of  r.  ,  and  (j) .  ,  . 

The  initial  values  of  tp   are  assigned  much  more  simply. 
Since  we  already  know  that 

^  =  H  , 

is  a  constant,  on  the  upper  boundary  and 

f   =   0 

on  the  lower  boundary,  we  approximate  f   by  the  linear  inter- 
polation 

^•,k  -  h  • 


8if 


4. 11  A  Description  of  the  Complete  Algorithm 

The  algorithm  proceeds  in  the  following  manner: 

1.  The  input  parameters  are  assigned. 

2.  The  initial  values  for  the  functions  i}/ ,    r  and  (j)  are 
computed. 

3.  The  time  is  assumed  to  move  forward  one  step.   The  values 
of  the  functions  Tp ,    r  and  \)   are  updated  at  all  the  interior 
points  of  R  according  to  equations  (34),  (56)  and  (37)' 

4.  The  values  of  r  and  ^   on  the  upper  boundary  of  R  are  up- 
dated by  iterating  equation  (42)  one  time.   This  applies  both 
the  tangential  shift  and  the  shift  to  keep  these  points  on  the 
fixed  boundary  of  D. 

5.  The  values  of  r  and  ^   on  the  lower  boundary  of  R  are  up- 
dated according  to  equation  (43).   This  consists  of  both  the 
tangential  shift  and  the  normal  shift. 

6.  The  approximate  solution  is  normalized,  according  to  the 
method  described  in  Section  7« 

7.  A  check  is  made  to  see  if  the  algorithm  has  converged. 
The  residuals  of  the  partial  differential  equations  (34),  (36) 
and  (37)  on  the  interior  of  R  as  well  as  those  of  the  boundary 
shifts  (42)  and  (43)  are  examined.   If  the  absolute  value  of 
each  one  is  less  than  some  arbitrary  pre-assigned  error  then 
convergence  has  been  obtained  and  we  stop  calculating;  otherwise 
we  go  back  to  step  3  and  repeat  the  process. 
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Chapter  5  -  Results 

5. 1   Introduction 

The  algorithm  described  in  Chapter  4  has  been  implemented 
on  the  CDC  66OO.   Various  geometric  configurations  for  the 
fixed  boundary  have  been  used  to  test  the  rate  of  convergence 
and  the  accuracy  of  this  method.   The  physical  parameters 
employed  were  chosen  to  be  near  those  of  the  experimental  machine 
SCYLLAC.   Our  main  result  is  that  with  appropriately  chosen 
weight  factors  the  method  converges  to  the  correct  solution  as 
the  artificial  time  t  goes  to  infinity  and  the  solution  of  the 
difference  equations  converges  with  second  order  accuracy  to  the 
solution  of  the  differential  equations  as  the  mesh  width  h  goes 
to  zero. 


5.2  Implementation  of  the  Method 

As  mentioned  above,  the  algorithm  was  written  in  FORTRAN 
and  run  on  the  CDC  66OO  of  the  AEC  Computing  and  Applied 
Mathematics  Center.   The  program  requires  about  250OO  decimal 
words  of  memory,  exclusive  of  the  storage  necessary  for  the 
functions  ijj{v.,y)  ,    r(u,v),  (t)(u,v)  and  c(u,v). 

As  it  is  presently  written,  the  program  is  not  as  efficient 
as  it  could  be.   In  order  to  decrease  programming  and  debugging 
time,  FORTRAN  was  used  as  a  programming  language.   This  choice 
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has  probably  increased  the  execution  time  by  a  factor  of  at 
least  two.   In  fact,  an  independent  study  by  R.  Hockney  [15]  has 
shown  that  an  Increase  in  efficiency  by  a  factor  of  five  or  six 
Is  possible  by  coding  in  machine  language  instead  of  FORTRAN  on 
the  CDC  6600. 

An  analysis  was  made  of  the  amount  of  time  the  computer 
required  for  each  segment  of  the  program.   The  results  are  given 
in  Table  I . 

In  addition,  a  version  of  the  program  was  written  for  use 
with  an  experimental  time-sharing  system  which  was  available. 
This  version  of  the  program  proved  to  be  of  aid  during  the  period 
when  the  proper  values  of  the  weight  functions  were  being 
selected. 


5. 3  Selection  of  the  Physical  Parameters 

The  physical  parameters  which  were  chosen  to  be  used  in 
the  actual  computer  runs  were  based  on  those  of  the  device 
SCYLLAC .   First,  the  twist  k  was  taken  to  be  one,  as  was  Be,,  the 
component  in  the  ignorable  direction  of  the  magnetic  field  B. 
Thus,  the  partial  differential  for  ^  was  brought  into  the  scaled 
form  of  equation  (I.50).   Second,  it  was  found  that  the  deter- 
mining factor  for  the  location  of  the  free  boundary  is  the  ratio 
^   of  the  periodic  strip  R  in  the  (u,v)  plane.   This  was  to  be 
expected,  since  ^^  is  the  Riemann  modulus  of  the  conformal  mapping 
It  was  chosen  so  that  the  ratio  of  the  radius  of  the  fixed  bound- 
ary to  that  of  the  free  boundary  would  be  approximately  five  to 
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one.   Again,  our  selection  was  based  on  the  parameters  of  the 
machine  SCYLLAC.   Finally,  the  parameter  H,  the  value  of  the 
flux  function  on  the  fixed  boundary,  was  arbitrarily  defined  to 
be  one. 


5.^  A  Description  of  the  Cases  Run 

The  first  case  run  was  that  of  the  magnetic  field  with 
circular  outer  boundary  defined  by 

F(r,(t))  =  r-  1  =  0  . 


Our  aim  in  this  case  was  to  check  the  accuracy  of  the  computed 
solution,  since  the  exact  solution  is  already  known  (cf.  Section 
1.9)'   The  program  was  run  with  the  mesh  sizes  h,  p-,  -jj-  and  q- > 
where  h  =  A,    until  all  the  residual  errors  were  less  than  10"  . 
It  was  verified  that 

1.  The  position  of  the  free  boundary  approached  the  exact 
solution  with  truncation  error  of  the  form  0{h  )  as  the  mesh 
width  h  went  to  zero. 

2.  The  complex  function  f(z),  which  is  theoretically  zero, 
did  approach  zero  with  error  0{h  ). 

3.  The  inverse  Beltrami  equations 


and 


2 

were  satisfied  to  order  h  . 

H- .      When  the  relaxation  numbers  a.  were  chosen  according  to 
formula  (4.46),  the  number  of  cycles  necessary  to  achieve  con- 
vergence doubled  as  the  mesh  width  was  halved.   Therefore,  the 
rate  of  convergence  is  O(r-)  and  the  linear  theory  of  Section 
4.8  is  applicable. 

Table  II  summarizes  the  results  of  this  case. 
The  second  case  run  was  the  ellipse 

F(r,(t))  =  (r  cos  (j))  +p  (r  sin  (j))  -1  -  0  , 

where  p  =  .75-   While  this  case  is  still  fairly  simple,  the  solu- 
tion cannot  be  obtained  explicitly  and  the  variables  do  not 

2 

separate.   Again,  the  accuracy  was  observed  to  be  of  order  h  . 

Table  III  and  Figure  5-1  summarize  the  results  of  our  second 
calculation. 

The  third  case  which  we  ran  had  an  outer  boundary  specified 

by 

(1)     F(r,(|))  =  (r  cos  (j))  +(r  sin  (j))  -1  =  0  . 

Clearly,  equation  (l)  is  equivalent  to  the  formula 

X  +y  =  1 

in  rectangular  coordinates.   The  interesting  aspect  of  case 
number  three  is  that  at  the  points 

1    _   TT      57r     57r     Jv 
4>   -  If  ^  -IT  '  T" '  T" 
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the  outer  boundary  curves  sharply,  behaving  approximately  as 
if  at  a  corner.   Since  the  mapping  is  conformal,  the  images  of 
the  mesh  points  on  the  upper  boundary  of  R  tended  to  cluster 
away  from  these  four  boundary  points.   However,  when  the  total 
number  of  mesh  points  exceeded  2000,  the  convergence  and 
accuracy  obtained  were  satisfactory  (cf.  Figure  5' 2). 
The  fourth  case  run  was  the  figure 

P  2     2 

F(r,(l))  =  (r  cos  (j)  -  .1)  +  (r  sin  ^)     -   1.5 

This  configuration  is  of  great  practical  interest,  since  it  is 
one  of  those  being  considered  for  the  device  SCYLLAC  which  is 
currently  under  construction  at  Los  Alamos.   It  was  observed 
that  the  shape  of  the  free  boundary  followed  that  of  the  fixed 
boundary  (cf.  Figure  5'3)' 

The  last  case  which  we  ran  had  the  fixed  boundary  specified 
by  the  equation 

F(r,(t))  =  (r  cos  ^-1,5)  +  (r  sin  \>)     -1  =  0, 

which  describes  a  circle  whose  center  Is  displaced  1.5  units 
from  the  origin  of  the  coordinate  system  (cf.  Figure  5'^)-   When 
examined  in  three  dimensions,  the  magnetic  field  produced  is  a 
helix  which  spirals  around  the  z-axls.   The  field  has  the 
property  that  the  singularity  of  the  partial  differential  equa- 
tion that  occurs  at  the  origin  of  the  coordinate  system  is  now 
located  outside  the  fixed  boundary,  whereas  In  all  previous 
cases  it  was  located  within  the  plasma. 
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From  a  mathematical  point  of  view,  this  case  proved  to  be 
the  most  difficult.   The  (r,(j))  coordinate  system  is  not  well 
suited  to  the  geometry.   Moreover,  the  images  of  the  mesh  points 
in  the  (u,v)  plane  tended  to  cluster  at  the  places  where  the 
two  boundaries  are  closest.   Further  refinement  of  the  mesh 
would  not  be  very  practical,  since  the  cost  would  be  high  and 
the  answer  has  already  been  obtained  with  sufficient  accuracy 
over  most  of  the  domain  D.   A  better  method  would  be  to  sub- 
divide the  auxiliary  domain  R  into  regions,  with  each  region 
having  its  own  mesh  size.   Thus,  a  more  refined  mesh  could  be 
imposed  on  any  region  of  R  which  maps  onto  a  larger  section  of  D, 
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Appendix  A  -  Conversion  Relations 


Cylindrical 
Coordinates 


Coordinates  for 
Helical  Symmetry 


basis  vectors: 

^r'  ^d'    ^z  ^r'  ^V  ^^ 

conversion  relations  between  basis  vectors: 
e  =  e^ 


r    r 

eg  =  e^+kre^ 

e^  =  -kre^+e^ 
Inner  products: 


e  ^e   =  1 

r  r 

e^*  en  =  1 


e  -e   =0 
z   z 


all  others  are  zero 


e   =  e 
r    r 


kre 


'i' 


\^  V      — 


o 

kren  +  e 
y   z 


e  •  e  =1 
r  r 

^   ^    o 
all  others  are  zero 


cross  products: 


^r^^e  =  ^z 


e  xe^  =  e, 


^^  =  ^C 


^e^^z  -  ^r 


e  X  e  =  e, 
z   r    ( 


e(i,^^C  =  ^ 


e»  X  e   =  e 


^ 


typical   vector; 


B  =  B  e    +  Boen  +B„e^ 
r  r        9   y        z   z 


B  =  Vr+ V't^^'^^^C 
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conversion  relations  between  components 


B  =  B 
r    r 

B I  +  krB^ 
B   =   'P   k 


B   =  B 
r    r 


Bi  =  Bfl  -  krB 
(f)     9      z 


B   = 
z 


krBi  +  B^ 


a 


B^  =  krBa  +B 
C      9    z 


operator  V: 


-r  Br    r  M   ®z  ^ 


^  =  ^r  ^  +  ^  ^<t)  B^ 


gradient; 


r  r  r    r   0    z^z 


Vtl/  =  e   ilJ     +  -^^  e 


divergence : 

r        Sr  r  69        dz 


1    ^(^B    )  SB, 

V.B   =   -  ^-^—  +  -  ^ 

r        or  r  d(p 


curl : 
VxB 


/•■      BB        hB, 

el-       ^  - 
r\  r   69        dz 


+   e 


SB        SB 
r  _    z 

9l^'5i~    BF" 


+  e 


z   r 

< 


S(rBg) 
Sr 


,    SB 
1        r 


Vx  B 


e   fl 


SBp\  r        2    SB 

2    :.    ,rBi\      2kBJ' 

+    r    ^'v^^;  ^      2 
a  a 
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Appendix  B  -  Tables  and  Figures 


Table  I 


Where  the  program  spent  its  time 


Section  of  program 


Percentage  of  time  spent 


Calculation  of 
coefficients 

Iteration  of  equation 
for  psi 

Iteration  of  equation 
for  R 

Iteration  of  equation 
for  theta 

Calculation  of  the 
shift  for  the 
free  boundary 

Calculation  of  the 
shift  for  the 
fixed  boundary- 
Normalization  of 
the  solution 

Library  routines, 
such  as  square  root 

Other  routines 


28,1 
26.9 

17.3 
8.2 

1.7 


8.9 
5.3 

3.0 

100.0 


9^ 


Table  II 


Circular  Fixed  Boundary 


Basic  mesh  width:   h 


Convergence  parameters 


Radius  of  the 

fixed  boundary: 

r  = 

1 

0 

Auxiliary  domain: 

Period 

A  = 

6 

8 

Height 

B  = 

2 

0 

a^  =  3.0 

P,  =  .6 

a^  =  8.0 

^2=  -^ 

a^  =  2.7 

p^  =  .6 

P4  -  -6 

mesh 
width 

h 

h/2 

h/4 

h/8 

number  of 

interior 

points 

68 

306 

1292 

5304 

number  of 
iterations 
for 
convergence 

392 

459 

699 

1120 

radius  of 
the  free 
boundary 

.165956 

.190293 

.194380 

.195324 

value  of  the 

2 
parameter  A 

4,098699 

3.262929 



3.160384 

3.137477 
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Table  III 


Elliptical  Fixed  Boundary 


Basic  mesh  width: 

h 

= 

A 

Convergence 

parameters: 

Maximum  error: 

E 

= 

10-5 

a^  ==  3.0 

^-^  =   -6 

Auxiliary  domain: 

Qp  =  8.0 

^2=  .2 

Period 

A 

= 

6.8 

a^  =  3.5 

p^  =  .6 

Height 

B 

= 

2.0 

Pi+  -  .6 

mesh 
width 

h 

h/2 

h/4 

h/8 

i 

number  of 

Interior 

points 

68 

306 

1292 

5304 

number  of 
iterations 
for 
convergence 

312 

382 

554 

914 

i 

radius  of  the 
free  boundary 
measured  along 
the  X-axis 

.188449 

.210185 

.214308 

.215318 

value  of  the 
parameter  A 

2.430560 

2.110442 

2.065419 

2.054323 
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FIGURE  5.1  -      ELLIPTICAL  OUTER  BOUNDARY 
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FIGURE  5.2 
FIXED  BOUNDARY  WITH  ROUNDED  CORNERS 
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FIGURE  5.3  -  SLIGHTLY    OFFSET  OUTER  BOUNDARY 
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FIGURE  5.4    HIGHLY  OFF-CENTER  OUTER  BOUNDARY 
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*0F 

♦  OF 
♦IF 
♦IF 
♦Bl 
COM 
REW 
#HE 
REQ 
UPD 
RET 
♦TI 
REW 
RUN 
REW 
COP 

♦  GO 

♦  CX 
ERE 

♦  TE 

♦  A 
REW 
MAP 
♦IF 
COM 
REW 
COP 
RET 
LGO 
EXI 

♦  FI 
♦TI 
REW 
REW 
COP 
REW 
COP 
♦FI 

♦  IF 
COM 
REW 
RUN 
GRP 


0107,T177,CM125000,   FRIEDMAN  N, 

FSWl, 

FSW2, 

(SW1)/*GOTO/*R1/ 

< , NOT, COMMON  TAPllO 7 ) REW I ND( TAP110  7 ) 


MON 

IND( 

G 

UEST 

ATE( 

URN( 

ME. 

INDC 

(S,« 

IND< 

YCF( 

TO/^ 

IT. 

DIT, 

RM. 

IND( 

(ON) 

(  ,N0 

MON 

IND( 

YBF( 

URN( 

(TAP 

T. 

N. 

ME. 

IND( 

IND( 

YBF( 

IND( 

YBF( 

N. 

(  ,N0 

MON 

IND< 

(S,, 

1107 


TAP1107, 
TAPH07) 

DRUMsNEWPL, 
N) 
NEWPL) 

COMPILE) 

, COMPILE, EROUT) 

EROUT) 

ERDUT) 

A/ 


LGO) 

T, COMMON    BARH07)REWIND(8AR1107) 

BAR1107, 

8AR1107) 

LGO,BAR1107) 

BAR1107) 

1107) 


TAPE65) 
TAPE90) 
TAFE90> 
TAPE89) 
TAPE89) 

T, COMMON  QRP1107>RewlND(QRPli07) 

GRP1107, 

GRP1107) 

##|GRP1107) 


EOR 


•COMDECK,  LSTS 

COMMON/LISTS/HI,     H2,     H3,     HA,     HS 
UIMENSlOiNl    Hl(140>,    H2(140)#    H3(140), 


M4(140),     H    5(140) 
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DIMENSION  Z  (17) 

COMMON  /TWOPI/  TWOPI,  PlC 

COMMON/SLINF/ISWT,  Cl,  C2 

COMMON/ZET/ZETAi  ZETA2,  C 

COMMON/K/KAY,  KAY2,  FCON, 

REAL  KAY,  KAY2 

COMMON  /SWT/  IKONi  IBND, 

C0MM0N/ZARA/NC0N,NC0NIN,C 

REAL  NCON,  NCONIN 

INTEGER  TSFLAQ 


(140,41),  THETA  (140,41),  C  (140, 4i) 

,  R>,  (TH,  THETA),  (CH,  C> 
TH  (2),  CH  (2) 


ON 

,  ERATU,  ER4TV 
FFSET,  nFFS=T2 
KINA,  M2NA.  NA 

PN,  IPOINT,  JPOINT,  REXP 
OMN 


INTEGER  TSFLAQ 

COMMON  /TS/TSFLAG 
*C0MDECK,VARDIM 

DIMENSION  A  (140, 
♦DECK, TOR 

PROGRAM  T0R2 

1  (TAPE67S1002,  T 

2  TAPE89S402,  TAP 

C       SE 
C       SE 
C       SE 
C       Js  1  DOWN 
C       Ja  2   UP 
C       TAPE  65  IS 
C         TAPE  67 
C       TAPE89  IS 
C  TAPE  90  IS 

C  THE  MAIN  LOOP 

*CALL,CBLK 
DATA 


41) 


002,  INPUT8202,  OUTPUTalOO?,  TAPE65s402, 

,TFUT,  TAPE698nUTPUT) 

) INITIALIZE, 

,.  FILL  FRQvi  COMMON  FILE  TAPF67 
SAVE  PSI,  ?,  THETA,  ON  COMMON  FILE  TAPE6 


FOR  THE  RESIDUAL 

IS  FOR  PSI,  R,  Th 

FOR  R  (5)  AND  THE 

FOR  R  (5)  A 

IS  AT  ST 


S  OF  PSI,  R,  THETA 

ETA 

TA  (5)  ON  THE  FIXED  BOUNDARY. 

NC  THETA  (5)  ON  THE  FREE  BOUNDARY, 

ATEMENT  NUM3ER  5078 


DATA  NA,  MA/140,41/ 
COMMON/REL/WP,WH,WT 
COMMON  /SIZE/  M,  Ml 
C0MM0N/MESH/H2 
COMMON/CNT/IP, IT, IR 
COMMON/BND/IFX, IFR 
COMMON  /B/BTFX,  9F,  9TFR, 
COMMON  /TIME/  RP,  RK,  RT, 
COMMON/LAM/LAMBUA 
COMMON/DIMEN/NAl,  N 

or Al   I  AMRn A 


I  MPl,  N,  Nl,  NPl 


BNFR 
RFR,  RFX 


REAL  LAMBDA 
COMMON/CONST/CONST 
COMMON/AREA/AP,  B 
REAL  LDERB,  NDERB, 
INTEGER  RFLAG 
KOUIVALENCE  (Z 
EQUIVALENCE  (Z 


ASQ,  ^A 


(1), 
(10), 


LDERT, 
Zl) 


ZIO) 


NDERT 

(Z  (4),  Z4),  (  Z(7).  Z7) 

,  (Z  (12),  Z12),  (Z  (14),  Z14),  (Z  (16), 
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c 
c 


c 
c 
c 
c 
c 

R 
C 
C 

c 
c 

c 
c 


DIME 
DIME 
INTE 
REAL 
UATA 
UATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
UATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
I8ND 
IBiMD 
COMH 
DIHE 
I 
N 
I 

I 
P 

0 

c 


10/ 


.94/ 

1,0/ 


Z16) 

NSION  1ST  (8) 
NSION  inER  (36) 
GER  OLDN,  OlDM 

NEX,  M6X 

KAY,  ZETA/1, 0,1,0/ 

RAT/1.0/ 

Nl,  Ml/34, 

AP  /6,8/ 

GV,  QW/  ,87, 

CI,  C2/0,0,  _.  . 

CP,  CR,  CT/3,0,  3.3,  2,7/ 

Bl,  H2,  B3i  H4/.6, .2, .6, .6/ 

ERR/1. OE-8/ 

I*  III  12,  13,  14/7,50,  50,  50,  2500/ 

15/100/ 

JPOINT/0/ 

RP,  RR,  RTi  RFR,  RFX/  5  *  0.0/ 

IP,  IT,  IR»  IPX,  IFF/  5*0/ 

Rl,  R2,  R3/3HPSI,  ll-R,  5HTHETA/ 

OFFSET/0,0/ 

CON,  CONlN/1,0,  0.5/ 

IKOM,  IBND/0»1/ 

P/1.0/ 

IPOINT  /3/ 

TSFLAG/0/ 

REXP/n.O/ 

ISWT  /  3/ 
DO 
=  ii        BOauuiNi 
ON/FXCON/  JJ,  KK 
NSION  MSG  (3) 
NPUT  VARIABLES 
AMELIST  SWITCH 
KON  IF  ZERO  DO 

MAPPING 
BND   USE  A  FIXtD  BOUNDARY  DEFINED  BY 

(X  "  OFFSET)  **  IBND  ♦  (P*Y)  •♦  IBND  -  CON  *♦  iBND  =0 
FFSET 
ON 


1 

=  2 


BOUNDARY  1. 
BOUNDARY  2, 


FREE  EOLNDARY  PROBLEM,  IF  GT  ZERO  DO  A  COnFORMA 


IPOINT   SHOULD  BE  3 
ISWT    SHOULD  HE  3 
REXP   NOT  USED 


C 
C 
C 
C 
C 
C 
C 


NAMELIST  PARI 

KAY   THE  TWIST  K 

ZETA   THE  COMPONENT  OF  THE  MAGNETIC  FILED  IN  THE 

IQNORAflLE  DIRECTION, 
AP    THE  PERIOD  OF  THE  STRIP  IN  THE  (U,V)  PLANE, 
CI   THE  CONSTANT  VALUE  CF  PS  I  ON  THE  FREE  BOUNDARY. 
C2    THE  CONSTANT  VALLE  OF  PSI  ON  THE  FIXED  BOUNDARY 
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c 
c 
c 
c 


OFFSET,  CON   SEE  THE  CESCRIFTIOM  OF  THE  NAM6LIST  SWITCH, 

CONIN   USED  ONLY  IN  CCNFORMAU  MAPPING, 

Nl    NUMBER  OF  MESHPQINTS  ALONG  THE  U-AXIS, 

Ml   NUMBER  MINUS  ONE  CF  THE  MESHPOINTS  ALONG  THE  V-AXIS, 


C 

C 
C 
C 

c 
c 
c 
c 
c 
c 


NAM6L1ST  PAR2 

THE  WEIGHTS 

CP   C  FOR  PSI 

CR   C  FOR  R. 

CT   C  FOR  THETA, 

Bl   WEIGHT  FOR  TANGENTIAL  ShIFT  ON)  THE  FIXED  BnUMQARY, 

82   WEIGHT  FOR  THE  SHIFT  TO  KEEP  THE  FIXED  BOUNDARY 

IN  THE  CORRECT  POSITION, 
B3   WEIGHT  FOR  THE  TANGENTIAL  SHIpT  ON  THE  FMEE  BOUNDARY, 
84   WEIGHT  FOR  THE  NORMAL  SHIFT  OM  THE  FREE  BOUNDARY. 


C 
C 

c 
c 
c 
c 
c 
c 

c 
c 
c 
c 


7000 
7001 
7002 
7003 
7004 


170 


6666 


NAMELIST  PAR3 

I  SEEE  DESCRIPTION  AFTER  STATEMENT  190 
ERR  MAXIMUM  ERROR  FOR  CCNVERGENCE. 

II  AFTER  II  ITERATIONS,  PRINT  TH=  MAXIMUM  RESIDUALS. 
12,13    GOVERN  RATE  OF  AUXILLIARY  PRINTOUT, 

14  THE  MAXIMUM  NUMBER  OF  ITERATIONS  TO  BE  PERFORMAD  UNDER 
SCHEME  U. 

15  GOVERNS  RATE  AT  WHICH  THE  REMAlNNQ  TIME  IS  CHECKED, 


NAM 
NEX 
HEX 
B 
NAMELI 
NAMELI 
NAMELI 
NAMELI 
NAMELI 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
CALL  S 
CALL  S 
TLs  1 
TL»  ,0 
PRINT 
FORMAT 
CALL  T 
CONTIN 
PRINT 
FORMAT 
PRINT 
READ  S 
PRINT 


tLIST 
RATI 
RATI 
OTH  N 
ST  /S 
ST/PA 
ST/PA 
ST/PA 
ST/PA 
(*  R 
(♦  R 
(♦ 
<♦ 
(* 
TATIJS 
eCOND 
ST  (3 
01  ♦ 
170, 

(*  T 
SCHCK 
UE 
6666 

(1H3 
7000 
WITCH 
6667 


t^a/  1  ,  cnn  ,  1  J. ,  i  ^  ,  1 

R4/NEX,MEX,QV,Qk< 

EAD  SWITCH*) 

PAD  PARI*) 

EAD  PAR2*) 

EAD  PAR3*) 

EAD  PAR4*) 

(1ST) 

(TX) 
) 

TL 
TL 
IME  REMAINING  IS*, 


•  10,0  ♦  TX 


F8,3,  ♦  SECONDS,*) 


/IHl  ) 
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6667  rO«hAT  <*  I  A^  T0R2 *) 

PRINT  6668,  IKON,  IBND,  F,  IPOLMT*  R=XP 

6668  FORMAT  (IH  ,  ♦IKONs*,  13,  *ieNr=*#  H,  *Pa*,  VH,A, 

1   *iPOiNT=*,  13,  *REXPs*,  re, 4) 

PRINT  7001 
READ  (5,  PARI) 
PRINT  700? 
READ  (5,  PAR2) 

0rFSET2s  OFFStT  **  2 

TWOPIs  2,0  ♦  3.14159265358979 

PICONS  TWOPl 

P=  ABS  (P) 

PNs  P  ♦♦  I8ND 

ZETA2=  ZETA  ♦*  2 

KAY2=  KAY  *♦  2 

FCONs  -  2.0  ♦  KAY  ♦  ZETA 

CONSTs  CON 

NCONs  CON  *♦  IbNU 

NCOMNs  CONIN  ♦* 


60 
61 

4657 


IBND 

IF  (  IKON  ,NE.  0>PKINT  60 
FORPAT  C*  DO  CONFORHAL  MAPPING.*) 
PRINT  61,  ZETA,  KAY 
FORMAT  (*  PARAMETERS, */10X»  *ZETAs*, ri2, 6/lOX,  *Ka*,  F12.6) 


104 
101 


SIZE 

0 
RECTANGLE 


IS*,  F10,7  ) 


IS*,  Fin, 6,  *  BY*,  F10.6  ) 
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CONST  IS  THE  VALUE  OF  R  AT  Vai, 

PRINT  4657,  CONST 

FORMAT  (*  THE  INITIAL  FIXEC  GOLNDARY  IS  AT**  FiO.6) 
SET  UP  SIZE  OF  THE  ARRAYS, 

H=  AP  /  Nl 

H2=  2.0  *  H 

a=  Ml  *  H 

PRINT  104,  M 

FORMAT  (*  MESH 

PRINT  101,  AP, 

FORMAT  (  *  THE 

N=  Nl  ♦  1 

NPls  N  ♦  1 

IF  (NPl  ,GT.  NA>  QO  TO  401 

Ms  Ml  ♦  1 

MPls  M  ♦  1 

NAla  NA  -  1 

NASQs  Ml  *  NA  ♦  N 

MlNAa  Ml  *  NA 

M2NAS  MINA  -  NA 

IF  (M  ,GT,  MA)GO  TO  402 

NEQs  (Ml  -  1)  «  Nl 

FORMAT  (*  Ca*.  3F10.5) 

Ts  SORT  (AP  *  B> 

WTa  RELAX  (CT»  H  /  T) 

WRa  RELAX  (CR,  H  /  T) 

WPa  RELAX  (CP,  H  /  T) 
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110   FORMAT  (*  RELAXATION  FOR  PS  I  IS*,  F10,6,  •   FOR  R  IS  *,  FIO.61 
1    *   FOR  THETA  IS  *,  FlO.6) 
BTFX*  H  *  Bl 
BF-  H  ♦  B2 
BTFRa  H  *  B3 
8NFR3  H  *  B4 
IF  (SENSE  SWITCH  1)820,819 

819  CAUL  INIT 
PRINT  830 

630   FORMAT  <♦  INITIALIZE  ALL  VARIABLES.*) 
GO  TO  840 

820  CONTINUE 
PRINT  835 

835   FORMAT  (♦  FILL  VARIABLES  FROM  COMMON  FILE  TAPE67,*) 
IF  (OFFSET  .GT.  CONST)PICONa  0,0 
CALL  FILL 

821  CONTINUE 
840   CONTINUE 

DO  87  Ja  1,  NPl 

00  87  Ks  1,  M 

C  (J,K)a  CALCC  (J,  K) 
87    CONTINUE 
90    CONTINUE 

PRINT  1201 

1201  FORMAT  (*  WHAT  KIND  OF  A  PROBLEM  IS  THIS,*) 
IF  (IKON  ,En,  0)GO  TO  1159 

PRINT  1203 
GO  TO  1200 

1199  PRINT  1202 

1202  FORMAT  (lOX,  *A  FREE  BOUNDARY  PROBLE^l,*) 

1203  FORMAT  (lOX,  ♦A  CONFORMAL  NAPPING,*) 

1200  PRINT  1204 

1204  FORMAT  (*  WHAT  IS  THE  MESH  SIZE.*) 
PRINT  1302,  H 

1302  FORMAT  (lOX,  ♦Ha*,  F12,8) 

PRINT  1301 
1301  FORMAT  (*  WHAT  IS  THE  GRID  SIZE,*) 

PRINT  1205,  Ml,  Nl 

1205  FORMAT  (lOX,  13,  *  BY*,  13) 
PRINT  1206,  NEQ 

1206  FORMAT  (*  HOW  MANY  EQUATIONS  ARE  THE^E.*/  lOX,  17) 
PRINT  1207 

1207  FORMAT  (*  WHAT  ARE  THE  PARAMETERS,*) 
PRINT  1208,  KAY,  ZETA 

1208  FORMAT  (lOX,  *Ka*,  FlO,4/  lOX,  *ZETAi*,  FlO,4) 
PRINT  1209 

1209  FORMAT  (*  WHAT  ARE  THE  IMTlAL  ROUiMDaRIES  OF    PSI.*) 
PRINT  1210,  C2,  CI 

1210  FORMAT  <inx,  *PSI  ON  THE  OLTER  BOUiJOARY  IS*,  Fin, 4/ 
1    40X,  *PSI  ON  THE  INNER  60LNCARY  IS*,  FlO.4) 

PRINT  1220 
1220  FORMAT  (*  WHAT  IS  THE  INITIAL  FIXED  30UNDARY,*) 
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1221 

1222 

1223 

1211 

1212 
1 

112 

111 

1 
2 
3 

WRITE 

690 

1 

692 


695 


696 


190 
205 


PRIM 

FORMAT 

PRINT 

FORMAT 

PRINT 

FORMAT 

PRINT 

FORMAT 

PRINT 

FORMAT 

lOX 
PRINT 
FORMAT 
PRINT 
FORMAT 

lOX 

♦  F 
lOX 

OUTPU 
WRITE 
FORMAT 

♦  6 
WRITE 
FORMAT 
ENDFIL 
PRINT 
FORMAT 

♦  N 
* 
*FR 

PRINT 
FORMAT 
4{* 
JJs  0 
KKa  LL 
TMs  TI 
NNrO 

PRINT 
READ  ( 
PRINT 
FORMAT 
IF  (I 
GO  TO 

512 
1=1   C 


1221, 
(lOX 

1222 
(*  W 

1223, 

(inx 

1211 

(*  w 
1212, 

(lOX 
,  *CP 
112 

{*  w 

111. 

(*  F 

,  *SH 
REE  B 
,  -NU 
T  HEA 
OUTPU 

(*1I 
X,  *R 
OUTPU 

(/// 
E  68 
695 

(IHl 

PSI 

FIX 

EE   N 

696 

<1H 


R  (1,  M) 
,  *R'»*,  F20,10  ) 

HAT  IS  THE  OFFSET,*) 

OFFSET 
,  ♦OFFSETS*,  F20.10  ) 

HAT  ARE  THE  RELAXATIQ;\)  ^JlJ^HPRS,♦) 

WP,  WR,  WT,  CP,  CR,  CT 
,  *WP3*,  F10,4,  *  WRa*,  FlO,4,  ♦  WTa*,  F10.4/ 
«♦,  FIG. 4,  *  CR=*,  F10,4,  *  CTs*.  FlO,4) 

HAT  ARE  THE  WEIGHT  FACTORS,*) 

rfTFX,  bF,  HTFR,  MNFR 

IXED  BOUNDARY,*/  lOX,  *TA\)  SHIFT  FOR  IS*.  FlO.6/ 

IFT  FOR  -F  IS*,   F10.6/ 

nuNDARY,*/   lOX,  *TAN  SHIFT  IS*,   F10,6/ 

RMAL  SHIFT  IS*,   FID, 6) 

DINGS 

T  TAPE  9(),  690 

TERATIONS*,10X,  *R  (5)*,  lOX,  *THETA  (5)* 

DOT  (5)*,  7X,  *THETADOT  (5)*,  3X,  *LAMqDA*) 

T  TAPE  89,  692 

*  FIXED  BOUNDARY.*) 


RESIDUAL  AT 
ED  TAN  SHIFT 
ORM  SHIFT*   ) 


R  RESIOUAL   AT       THETA  RES,   AT*. 
-F  FREE   TAN  SHIFT   *, 


,  3(*- 


2X,  2^--) 


II-  )) 


1=2 
1  =  3 
1  =  4 
1  =  5 
1=6 
1=7 
1  =  6 


a  MMs  0 
=  0.0 

7003 

5,  PANS) 

205,  I,  ERR 

(  //  *  DO  ROUTINE*, 

.EQ.  0)  CALL  EXIT 

(501,  502,  503,  504, 

,  513,  514,  515,  516, 

0  NOTHING 

0  NOTHING 

0  NOTHING 

0  NOTHING 

0  NOTHING 

RINT  OUT  THE  CURRENT  VALUES  OF 

TERATE  ON  THfc  INTERIOR  AND  THE 

TERATE  ON  THfc  INTERIOR  ONLY, 


14,    *    eRRs#,    El?, 5) 

505,    506,    507,    508, 
517,    513,    519,    520, 


^^0  9, 
5?1 


510, 
522, 


511, 
523) 


PSI,    R,    TrlETA, 
BOUNDARY. 


AND    C, 
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501 
502 
503 
504 
505 
506 


507 
5079 


1*9 

IslO 

1 

2 

3 

4 

5 

Isll 

1=12 

I3l3 

1  =  14 

1=15 

1  =  16 

1=17 

Isia 

1  =  19 

I  =  iJO 

1=21 

1=22 

1=23 

GO  TO 


I.E., 

ITERA 
READ 
THE  F 
DO  PS 


DO 
DO 
DO 
DO 


R 

TH 

FI 

FR 


TO 
TO 
TO 
TO 


206 


207 


1 
5070 

FO 


GO 

GO 

GO 

GO 

CALL 

CALL 

CALL 

GO  TO 

IBRb 

CONTI 

I5lB 

PRINT 

FORMA 

/♦ 

/* 

/♦ 

/* 

/♦ 

PRINT 

FORMA 

) 

ILIM3 

E=  0, 
GO  TO 
R  I  R 

1=1 
1  =  2  F 
1  =  3  D 
1  =  4  D 


READ 
ASSIG 
REFIN 
RESTA 
SAVE 
READ 
OBSOL 
C8S0L 
OBSOL 
PLOT 
DO  A 
MULTI 
PRINT 
190 
190 
190 
190 
190 
PUTOU 
NICEP 
PUTOU 

190 
1 

NUE 
(15  ♦ 
206, 
T  (* 
PRIN 
PLOT 
THE 
CHEC 
THE 
207, 
T  (* 


DO  A 
TE  ON 
AN  IT 
ORMAT 
I  EQU 
EOUAT 
ETA  E 
X6D  B 
EE  BO 
NEW  W 
NAN 
E  THE 
RT  TH 
THE  C 
VALUE 
ETE, 
ETE, 
ETE. 
A  QRA 
LINEA 
PLY  A 

OUT 


DIRIC 

THE  B 

ERATIO 

IS... 

ATION. 
ION. 

QUATIO 
OUNDAR 
UNDARY 
EIGHT 

EW  VAL 
MESH, 
E  PROG 
URRENT 
S  OF  P 


HLRT  PROBLEM, 

OUNCARY  ONLY, 

N  SEQUENCE  IN  Qvj  THE  NEXT  INPUT  CARD, 

BLANK,  NUMBER,  BLANK,  NUMBER, ETC. 


N. 
Y. 

FACTORS. 

UE  FOR  THE  VALUE  OF  R  ON  THE  FREE  dOUNUARY. 

RAM, 

VALUES  OF  PSI,9,THETA  ON  TAPE. 
SI,  R,  THETA  FROM  TAPE, 


PH  OF  THE  BOUNDARIES  ON  THE  LINE  PRINTER. 

R  TRANSFORMATION  OF  T-lE  PROBLEM  ALONG  THE  X-AXlS. 

LL    OF    PSI    BY    A    CONSTAnjT. 

CAUCHY  RIEMANN  CONDITIONS  ON  BOTH  ROUNDARIESi 


T  (PSI,  3HPSI) 

UT 

T  (C,  IHC  ) 


70)  /  (Nl  *  (M  -  1)) 

Il#  12.  13,  14,  151,  IPOIvjT  ,  JPOINT 
PRINT  RESIDUALS  EVERY*.  15,  *  TIMES,* 
T  R  (5)  AND  THETA  (5)  EVERY*,  15,  *  TIMES,* 

ERROR  GRAPH  EVERY*,  15,  ♦  TIMES,* 
MAXIMUM  NUMBER  OF  ITERATIQviS  IS*,  15, 
K  THE  TIME  EVERY*,  13,  *  ITERATIONS,* 
NORMALIZING  SChEME  IS  NUMBER*,  12  ,  *  AND*,  12  //) 

ISWT 
ISWTs*,  13,  *  THIS  NLMaER  SHOULO  BE  POSITIVE,  BEST  IS  3* 


14  ♦  JJ 
0 

(5078,  5100,  5082,  5092),  IBR 
EAD  IBR, 
DO  ALL, 

OLLOW  LIST  IDEA, 
0  ONLY  INTERIOR, 
0  ONLY  BOUNDARY, 
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I8R=5      ADJUST 
5078    COiMTINUE 

CALL  PSINT 

CALL  RINT 

CALL  THEINT 

CALL  FIXBMD 

CALL  FREBND 


FIXED  BOUNDARY 


5077 


CONTINUE 
Ea  AMAXl 


(  ABS  (Z  (1)) 


1  ABS 
KK=  KK 
IF  (KK 
JJa  JJ 
KKa  0 

CALL  NONLIN 
OMEGAS  2.0 


(Z  (12))# 

*  1 
.LT. 

♦  KK 


ABS  (2 


ABS  (Z 
(X4)), 


(4) ), ABS  <Z  (7)  ), 
AaS  (Z  a6))> 


ABS  (Z  (10))i 


ID  QO  TO  301 


(ALPHA) 

/  (1,0  ♦  ALPHA) 
1)5001,^002,5001 
JJ,  Z,  QHEQA 

,  M,  3(E10.1, 
12,  IH) 


2H  (, 

JJ.  Z, 
,  14, 
IH)  ), 


OMEGA 
3(E11,1, 
4(E11,1, 


.LE,  CONST)  CALL  OVERLAP 
12)  GO  TO  302 


( JJ) 


IF(TSFLAG  - 

5001  PRINT  7111, 
7111  FORMAT  (IH 

1    4(E10.1, 
GO  TO  5003 

5002  PRINT  711, 
711   FORMAT  (IH 

1    IH,  ,12, 

5003  COMTINUE 
IF  (OFFSET 

301  LLa  LL  ♦  1 
IF  (LL  .LT. 
LLs  0 
CALL  FRQUT 

302  MMs  MM  ♦  1 
IF  (MM  ,LT 
MMa  0 
XJJb  JJ 
WRITE  TAPE65, 

5071  NNr  NN  ♦  1 
IF  (NN  .LT.  I51)G0 
CALL  STATUS  ( 1ST) 
IF(  IST(3)  .LT,  10000) 
NNaO 

5072  IF  (JJ  ,GE 
IF  (E  ,GT, 
IF  (E  ,GT, 
JJs  JJ  ♦  KK 

KKs  LL=  MMs  NNs  0 
CALL  NONLIN  (ALPHA) 
OMEGAS  2,0  /  (1,0  ♦  ALPHA) 
PRINT  711,  JJ,  Z,  OMEGA 
CALL  FROUT  (JJ) 


2H 
), 


(  ,  12 
F9,4) 


1H< 


12,  IH)  ),/ 


2H 


( 

2H 


(, 


12, 

12 


IH)  ),  F'9.4  ) 


13)  GO  TO  5071 


XJJ,  Zl,  Z4,  Z7,  ZIO,  Z14,  Z16 


TO  5072 


GC  TO  403 


ILIM)GO  TO  190 
1.0E4)  GO  TO  404 
FRR)  GO  TO  5070 


904   FORMAT  (  *  TIME 
CALL  DONE  (JJ) 

Jj»n 


IS  ♦,  5F15,a  ) 
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610 

6100 

6101 

508 

5081 

5082 

509 

5091 

5092 

510 
9510 
715 
717 

5100 

5101 

5102 

5103 

5104 

5105 

5108 

5109 

511 

211 

511U 


GO  TO 

THE  R 
WRb  QV 
PRINT 
FOR^^AT 
PRINT 
FORMAT 
QO  TO 
IBRb  3 
DO  508 

Z  ( 
QO  TO 
CALL  P 
CALL  R 
CALL  T 
GO  TO 
IBRs  4 
DO  509 

Z  ( 
QO  TO 
CALL  F 
CALL  F 
GO  TO 
IBRb  2 
DO  951 

Z  ( 
READ  7 
FORMAT 
PRINT 
FORMAT 
IF  (ID 
QO  TO 
DO  510 

JP  = 

IF 

60 
CALL  P 

GO 
CALL  R 

GO 
CALL  T 

GO 
CALL  F 

GO 
CALL  F 

CON 
QO  TO 
PRINT 
FORMAT 
READ  ( 
RFLAG= 
PRINT 


190 

ELAXATION  FOR  R  IS  TCO  HIGH, 

♦  OMEGA 
6100 

(*  INCREASE  C«SUB«R.*) 
6101,  WR 


<*  OMEGAS 
5070 


•»  FID. 4) 


1  JXa  10,  17 

JX)s  0,0 

5079 

SINT 

INT 

HEINT 

5077 

1  JXs  1,  9 
JX}3  0,0 
5079 

IXBND 
REBND 
5077 


0  JXs 
JX>s 
15,  I 

(361 
717, 

(♦  I 
ER  (1 
5079 
8  Js 

IDER 
(JP  . 
TO  (5 
SINT 
TO  51 
INT 
TO  51 
HEINT 
TO  51 
IXBND 
TO  51 
REBND 
TINUE 
5077 
211 

(*  I 
5,  PA 
0 
103, 


I,  17 

0,0 
DER 
2) 

IDER 

DER  IS*,  4UI2) 

)  ,Ei3,  0)  GO  TO  190 

1,  36 
(J) 

EQ,  0)  GO  TO  5109 

101,  5102,  5103,  5104,  5105)   ,JP 


08 

ne 

08 
08 


NPUT  NEW  RELAXATION  FACTORS  VIA  PAR2*) 
R2> 

rp,  CR,  CT 


iin 


WT 


Ts  H  /  SQHT  (AP  *  B) 
WT*  RELAX  (CT»  T) 
WRs  RELAX  (CR,  T) 
WPs  RELAX  (rp,  T> 
PRINT  llOi  WP,  WR, 
BTrx=  H  ♦  Bl 
BF-    H  *  H2 
BTFps  H  ♦  B3 
BNFRs  H  *  B4 
PRINT  112 
PRINT  111,  HTFX, 
QO  TO  190 


BF,  BTFR,  RNFR 


SHIFT  FREE  BOUNDARY  IN  R, 

512  PRINT  212 

212   FORMAT  (*  SHIFT  FREE  BOUNDARY  IN  P,  INPUT  NEW  VALUE,*) 

READ  5120,  RX 
512Q  FORMAT  {F20,10) 

K»l 

DO  5122  Je  1,  NPl 
R  (J,  K)r  RX 
C  (J,  K>=  CALCC  (J,  K) 
5122     CONTINUE 

PRINT  680,  RX 
680   FORMAT  (*  THE  FREE  BOUNDARY  IS  NOW*,  FlO,6) 

GO  TO  190 

513  NEXs  MEXsl 
READ  PAR4 
OLUNs  N 
OLDMa  M 

Nls  Nl  ♦  NEX  ♦  ,5 
Mia  Ml  •  MEX  ♦  ,5 
Ns  M  ♦  1 
NPlB  N  ♦  1 
Ms  Ml  ♦  1 
MP1«  M  ♦  1 
NASQs  Ml  ♦  NA  ♦  N 
MlNAs  Ml  *  NA 
M2NAB  hlNA  -  NA 
NEQs  Nl  *  (  Ml  -  1) 
PRINT  1(16,  Ml,  Nl,  NEO 
106   FORMAT  (*  MESH  SUE  IS*,  316) 
IF  (NPl  ,GT.  NA)  GO  TO  401 


IF  (M 
Ms  AP 
H2»  2 
PRINT 
B-  Ml 
PRINT 
IF 
DO 


.GT, 
/  Nl 
0  *  H 
104,  H 
♦  H 
101, 
(OFFSET 
150  Jal, 


MA)  GO  TO  402 


AP, 
•  LT 


ti 


CON)GO  TC  151 


DO 


150  Ksl, 
X=  R  (J, 


OLDiM 
OLDM 
K)    «    COS 


(THETA    (J,K)) 


111 


c 
c 


150 
151 


Y  = 

H  < 

THE 

CON 

CONTIN 

CALL  S 

CALL  S 

CALL  S 

IF  (OF 


00 
DO 


175 
176 
FI 

19 
5137 


21 
22 


514 
515 


RE 


516 


517 

518 


519 

520 

521 


175 
175 
X  = 

Y  = 
P  ( 

THE 

CON 
CONTIN 
X  UP  T 
DO  19 

THE 
DO  22 

DO 


CON 
FIX  UP 
GO  TO 
GO  TO 
AVE  DA 
REWIND 
WRITE 
WRITE 
(((P 
END  FI 
PRINT 
GO  TO 
AD  DAT 
CALL  F 
PRINT 
GO  TO 
GO  TO 
CALL  N 
OMEGAS 
IF  (WR 
WRa  QV 
PRINT 
PRINT 
GO  TO 
GO  TO 
CALL  D 
GO  TO 
PRINT 


R  ( 
J,K 
TA 
TIN 
UE 
HIF 
HIF 
HIF 
FSE 
Js 
Ks 

R  < 
THE 
J«K 
TA 
TIN 
UE 
HET 
Ksl 
TA 
Jal 
21 
C 
C 
TIN 

PA 
511 
1 
TA 

67 
TAP 

TA 
SI 
LE 
903 
190 
A  F 
ILL 
835 
513 
190 

ONL 
2. 

.L 

• 

610 
610 
190 
190 
ONE 
190 
700 


J,K)  ♦  SIN  (THETA  (J,K)> 
)=  X 
<J,K)»  Y 

UE 


SI,    NEX,  HEX,  OLON,  O.DM) 

,      NEX,  KEX,  OLDN#  OuDM) 

HETA.  NEX,  fEX,  OLDN*  OlDM) 
T.  CON)  GO  TO  176 


J,    K) 

QRT  (X  **  2  ♦  Y  •* 

)«  ATAN2  <Y,  X) 


2) 


T  (P 
T  (R 
T  (T 
T  .L 
1,  N 
1,  M 
J,K) 
TA  ( 
)s  S 
(J«K 
UE 

A 

,M 

(NPl,  K>«  THETA  (2i  K)  -  PICON 

,  NPl 

Ka  1.  M 

(J,  K>s  CALCC  (J,  K) 
ONTINUE 
UE 

RAHETERS, 
0 

ON  COMMON  FILE  TAPE67, 

E  67,  N,  M 

PE  67, 

(J,  K),  R  (J,  K),  THETA  (J,  K)),  Ja  1,  NPl),  K*  1,  M) 

60 

1 

ROM  TAPE67 


IN  (ALPHA) 

0  /  (1.0  ♦  ALPHA) 

T,  QW  ♦  OMEGA)  GO  TO  190 

OMEGA 

0 

1,  WR 


(I) 


112 


522 

5222 

5221 
523 


401 
901 

402 
902 

403 
903 


4031 


ALU  DF  PSI  BY*.  riO.5) 


9031 
4032 
404 


4041 


READ  PA«1 
PRINT  PARI 

0FFSET2S  OFFSET  **  2 
IF  (OFFSET  .LT,  R  (1,1))  GC  TO  190 
PICONsO.O 
CALL  ORGSHF 
GO  TO  190 
CONTINUE 
READ  5120,  RAT 
PRINT  5222,  RAT 
FORMAT  (*  MULTIPLY 
DO  5221  Jal,  NPl 
UO  5221  Ksl,  M 

PSI  (J,K)=  RAT  ♦  PSI  (J,  K) 

CONTINUE 
GO  TO  190 
CALL  CRMANN 
GO  TO  190 

ERROR  MESSAGES  ARE  HERE, 
PRINT  901 

FORMAT  (  ♦  LENGTH  TOO  SMALL,*) 
CALL  EXIT 
PRINT  902 
FORMAT  (  * 
CALL  EXIT 
PRINT  903 
FORMAT  (* 
JJa  JJ  ♦  KK 

IF  (SENSE  SWITCH  2)4031,4032 
REWIND  67 

WRITE  TAPE  67,  N,  M 
WRITE  TAPE  67, 

(((PSI  (J,  K),  R 
END  FILE  67 
PRINT  9031 
FORMAT  (♦  SAVE  DATA 
CALL  EXIT 
PRINT  4041 
JJs  JJ  ♦  KK 
CALL  NONLIN  (ALPHA) 
OMEGAS  2,0  /  (1,0  ♦  ALPHA) 
PRINT  711,  JJ,  Z,  OMEGA 

FORMAT  (♦  THE  ERROR  IS  BECOMING  TOO  ,ARGE,*) 
GO  TO  190 
END 


WIDTH  TOO  SMALL.*) 


TIME  LIMIT, ,,.,,*) 


(Ji  K>,  THETA  (J.  K)),  Jsl#NPl  ),Kal,M) 


ON  COMMON  FILE  TAPE  67.*) 


C 
C 


SUBROUTINE  TSCHCK 

CHECK  TO  SEE  IF  IN 
SET  TSFLAG  TO  1  IF 


Tlh'E-SHARING    MODE, 

IN  BATCH  MODE,  OTHERWISE 


IN  TS  MODE 
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*CALL»CBLK 

CAUL  REQPP(4LRW6P  ♦  LOCF < TSFLAG ) ) 
IF  (TSFLAG  -  1)5,7,5 
7    PRIM  100 

100  FORMAT  (lOX,  ♦  I  AM  IN  BATCH  MODE,*) 
RETURN 

5     PRINT  101 

101  FORMAT  (lOX,  ♦  I  AM  IN  TS  MODE,*) 
RETURN 

END 


SUBROUTINE  NONLIN  (ALPHA) 

CHECK  THE  RELAXATION  FACTOR  FOR   ^, 
♦CALL,CBLK 

COMMON/SIZE/M,  Ml,  MPl,  N,  (Ml,  NPl 

AMAX=0,0 

DO  10  Ja2,  N 

DO  10  Kb2,  Ml 

T=  (ABS  (C  (J  ♦  1,  K)  •  C  (J  -  1,  K))  ♦  ABS  (C  <J,  K  ♦  1) 
1    -  C  (J,  K  •  1)))  /  C  (J,  K) 
IF  (T  .GT,  AMAX)  AMAX«T 
10       CONTINUE 

ALPHAS  AMAX  /  8,0 

RETURN 

END 


SUBROUTINE  INIT 

INITIALIZE  PSI,  THETA,  R, 
♦CALL.CBLK 

COMMON/SIZE/M,  Ml,  MPl,  h,    Nl,  NPl 

C0MM0N/MESH/H2 

COMMON/CONST/CONST 

COMMON/AREA/AP,  B 

Ha  H2  /  2,0 

Ss  C2  -  CI 

IF  (S)  18.  20,  22 

18  PRINT  19 

19  FORMAT  (*  THE  STREAMLINES  CF  PSI  ARE  IN  THE  WRONG  ORDER.  REVERSE 
ITHEM,*  ) 

TsCl 
CI"  C2 
C2«T 
GO  TO  21 

20  C1»0,0 
C2si,0 

21  Ss  C2  -  CI 
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22 
23 


FI 


25 


35 


40 
45 


55 


60 
65 


PRI 

FORM 

1    * 

LL  IN 

Ks  1 


DO 


4- 
00 


DO 


FILL 
DO 


DO 


ILL 
IF 


3 

301 

4 

500 

74 

75 


IF  ( 

Ks  M 

MTER 

HPs 

R  (1 
RXa 
DO  3 

P 

R  (1 

K=  K 
IF  ( 
IF  ( 

PRINj 
FORM 


QO 
DO 


NT  23,  C2.  CI 

AT  (*  THE  FiXfcD  BOUNCARY  IS  AT  PSIs*,  F8,4/ 

THE  FREE  BOUNDARY  IS  AT  PSIs*,  r8.4) 

PSI 

M 

5  Js  1,  NPl 
SI  (J,  l)s  CI 
SI  (J,  M)=  C2 

5  Ks  2,  Ml 

s  (K  -  1)  ♦  H 

SI  (J,  K)e  S  ♦  V  /  0  ♦  CI 

5  Ks  2,  Ml 

0  40  Js  2.  NPl 

PSI  (J,  K)»  PSI  (1,  K) 
ONTINUE 
IN  THETAs  -U 
Os  TWOPI  /  AP 

5  Js  1,  NPl 

=  (J  •  1)  •  H 

HETA  (J,  K)3  TWOPI  -  U  *  TO 

ONTINUE 

5  Js  1,  NPl 

0  60  Ks  2i  M 

HETA  (J,  K>s  THETA  (^,    X) 

ONTINUE 

N  Ra  EXP(V)  AND  C. 

IKON  ,GT.  0)QO  TO  101 

FILL  IN  R 

KAY  ,FQ,  0,0)  GO  TO  74 

1 

s  10 

H  /  MTER 
,  M>=  CONST 
R  <1#  K*l) 

01  KPsl,  MTER 

Xs  RX  •  HP  *  (TO  *  RX)  /  SCRT  (1.0  ♦  KAY2  *  RX  •*?) 
,K)s  RX 

-  1 
K  )  4  4  3 

R  (1,1)  .GT,  0,0)  GO  TC  110 
T  500 

AT  (*  ERROR  IN  INIT.  R  IS  NOT  POSATIVF,*) 
0  IIQ 
5  Ksl,  M 
3  (K-1)  ♦  H 
s  CONST  ♦  EXP  (-TO  ♦  8) 

(J,K)s  T  *  EXP  (TO  *  V) 
ONTINUE 
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101 


105 

110 


ao 

85 


GO  TO  110 

S-    CONST  -  CONIN 

DO  105  Kal,  M 

-  1)  ♦  H 
K)a  S  ♦  V 


Vs  (K 

R  (J, 

CONTINUE 

DO  85  K« 

DO  80 


/  8  ♦  CCNIN 


1,  M 

Js  2,  NPl 
R  <J,  K)*  R  (1,  K) 

CONTINUE 
CONTINUE 
IF  (CONST  .QT,  OFFSET)RETURN 

PICON*  0,0 
CAUL  ORQSHF 
RETURN 
END 


*CALL#CBl.K 

DO  300 
DO 


SUBROUTINE  ORGSHF 

TRANSLATE  THE  MAGNETIC  FIELD  AND 
COMMON  /SIZE/M,  Ml,  MPl,  N,  Nl,  NPl 


PLASMA  ALONG  THE  X-AXIS 


Ksl, 
Jb2, 
COS 
SIN 
R  (J, 


M 

N 
(THETA 
(THETA 
K) 


(J.K)) 
(J.  K)) 


300 
CT  = 
STs 
RT  = 
Ys  RT  ♦  ST 

XPs  RT  *  CT  ♦  OFFSET 
R  (J,  K)a  SORT  (XP  *♦  2 
THETA  (J,  K>B  ATAN2  (Y, 
300     CONTINUE 

DO  305  Kal,  M 

THETA  (1,  K)a  THETA  (N, 
R  (1,  K)a  R  (N,  K) 
THETA  (NPl,  K)a  THETA 
R  (NPl,  K)a  R  (2,  K) 
305      CONTINUE 

DO  400  Jsl,  NPl 

DO  400  Ksl,  M 
400      C  <J,K)a  CALCC  (J,K) 

RETURN 

END 


♦  Y 
XP) 


♦  •  2) 


K)  ♦  PICO^J 
(2,  K>  -  Pi:ON 


SUtJROUTINE  FILL 
FILL  PSI,  R, 
*CALL,CBLK 


THETA  FROM  MAGNETIC  TAPE 
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14 
15 
100 

200 


150 
390 


175 
395 

19 

20 


COMM 

COMH 

COMM 

REAL 

INTE 

REWI 

READ 

JT1» 

READ 

1   (( 

REWI 

TNb 

TMb 

NEXb 

MEXa 

IF  < 

IF  ( 

PRIN 

FORM 

PRIN 

FORM 

1    * 

GRXs 

IF  ( 

DO  1 

DO  1 

X 

Y 

P 

T 

C 

CONT 

CALL 

CALL 

CALL 

IF  ( 

DO  1 

DO  1 

X 

Y 

R 

T 

C 

CONT 

DO  1 

T 

RETU 

END 


QN/SI 
ON/ME 
QN/AR 
NEX, 
GER  0 
ND  67 
TAPE 
OLDN 
TAPE 
(PSI 
ND  67 
OLDN 
OLDM 
M  / 
Ml  / 
NEX  • 
HEX  - 
T  100 
AT  (* 
T  200 
AT  (* 
THE 
CFFS 
GRX  , 
50  Js 
50  Ks 
=  R  ( 
=  R  ( 
(  J.K 
HETA 
ONTIN 
INUE 
SHIF 
SHIF 
SHIF 
GRX  , 
75  J  = 
75  Ks 
=  R  ( 
=  THE 
(JiK 
HETA 
ONTIN 
INUE 
9  K  =  l 
HETA 
RN 


ZE/M,  Ml,  MPI,  N,  Nl,  MPl 

SH/H2 

EA/AP,  B 

MFX 
LDM,  OLDN 

67,  OLDN,  OLDM 
♦  1 
67, 
(J,  K),  R  (J,  K),  ThETA  (J,  K)),  J"  1.  JTl),  K=l,  OLDM) 

-  1 
'  1 

TN 

TM 

1.0)15,  14,  15 

1.0)15,20,15 

I  AM  CHANGING  THE  SIZE  OF  THE  INPUT  MATRICES,*) 
,  TM,  TN,  Ml,  Nl 

THE  OLD  SIZE  WAS*,  F4,0,  *  RY*,  F4,0/ 
NEW  SIZE  IS*,  14,  ♦  BY*,  14) 
ET  ♦*  IBND 

LT.  NOON)      GC  TO  3SC 
1,  OLDN 
1,  OLDM 
J,K)  •  COS 
J,K)  ♦  SIN 
)s  X 

(J,K)«  Y 
UE 


(THETA  (J,K)) 
(THETA  (J,K)) 


T  ( 
T  ( 
T  ( 
LT. 
1,  N 
1,  M 
J,K) 
TA  ( 
)-    S 
(J,K 
UE 


PSI,  NEX,  h'EX, 
R,  NEX,  h-EX, 
THETA, NEX,  KEX, 
NCON)      GC  Tn 


J.  K) 

QRT  (X  ♦♦  2  ♦  Y 

)s  ATAN2  (Y,  X> 


OLDN,  OlDM) 

OLD.vl,  O.DM) 

OLDN,  OuDM) 
395 


**  2) 


,  m 

(NPl,  K)s  THETA  (2,  K)  -  PI^ON 
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REAL  FUNCTION  RELAX  (C,  I- > 
RELAX*  2,0  /  (1,0  ♦  C  *  »•) 
RETURN 
END 


SUBROUTINE  SHIFT  (A,  MJ,  MK,CLCN,  OLDM) 
MOVE  ELEMENTS. 
•CALL.VARDIM 

COMMON/SIZE/M,  Ml,  MPl,  N, 
INTEGER  DELJ,  DELJ1#  DELK, 
INTEGER  OLDN,  OLDM 
REAL  MJ,  MK 
EXPAND 


11 


12 
13 


15 


16 
17 
21 

25 


Nl,  NPl 
DELKl 


MK  ♦  FLOAT  (K 


1) 


MJ 


*  FLOAT 
A{J,K> 


(J  -  I)  ♦  1,5 


<K 


KsOLDM 
KPa 
JsOLDN 
JP« 
A(JP,KP>s 
J»J-1 

IF(J)7,7,5 
KaK»l 

IF(K)9,9,3 
CONTINUE 
FILL  IN  HOLES, 
KHATaO 
UO  25  Kai,  OLDM 

KP=  MK  ♦  FLOAT 
JHATsO 
DO  21  J=l,  OLDN 

JPs  MJ  *  FLOAT  (J  -^  1) 
DELJa  JP  •  JHAT 
DELJla  DELJ  •  1 
IF  (DELJ1>13,13,11 
DO  12  LJ»  1.  DELvl 
A  (JHAT  ♦  Lw# 
A  (JP.KP))  /  DELJ 
CONTINUE 
DELKs  KP  -  KHAT 
DELKla  DELK  «  1 
IF  (nELKl)17,17,15 
DO  16  LJs  1.  DELv 
DO  16  LK=  1, 
A  (JHAT 
KHAT)  * 


1.5 


1)  ♦  1,5 


1.5 


KP)a  ((DELJ  -  LJ)  *  A  (JHAT,  KP)  ♦  LJ  * 


♦  LJ# 

JHATajP 
COiMTlNUE 
KHATaKP 
CONTINUE 


CELKl 
♦  LJ, 

LK  ♦  A  (JI-AT  * 

CONTINUE 


KHAT  ♦  LK)a  ( (DELK 
LJ»  <P>)  /  DELK 


LK)  *  A  (JHAT 
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FILL  IN  PEWIODIUITY 
DO  30  Ksl,M 
30     A(NP1,K)=A(2,K) 
RETURN 

END 


REAL  FUNCTION  CALCC  (J,  K) 
COHPUTE  THE  COEFFICIENT  C. 
*CALL,CHLK 

CALCCs  l.n  /  (R  (J,  K)  ♦*  2) 

RETURN 
END 


♦  KAY2 


REAL,  FUNCTION  OR 
CALCULATE  THE 
♦CALL,  CdLK 

CRs  -2,0  /  (R  (J, 

RETURN 
END 


(J,  K) 
DERIVATIVE 

K)  **  3) 


CF  THE  C:)EFFICIENT  C  XlTH  RESPECT  TO  R 


REAL  FUNCTION  CALCD  (J,  K) 
CALCULATE  THE  COEFFICIENT 
*CALL,CBLK 

Tl=  KAY2  ♦  H  (Ji  K) 

CALCDs  1,0  /  R  (J,  K)  •  (2.0 

RETURN 

END 


D, 


IX)    /  (1,0  ♦  Tl  *  R  (J,  K)) 


REAL  FUNCTION  CALCF  (J,  K) 

CALCULATE  THE  COEFFICIENT  F  ((J  IN  THE  PAPER). 
♦CALL.CBLK 

Tls  1,0  ♦  KAY2  *  (R  (J,  K)  *♦  2) 

CALCF=  FCON  /  Tl 

RETURN 

END 


SU8PQUTINE  FNC  (L»  F,  FR,  FThETA,  J,  K) 
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c 

C 

c 

L  =  0 

c 

LSI 

c 

U 

♦  CALL 

,CBLK 
EFFb 
IF  ( 
IF  ( 

501 

COST 
RX  = 
U-    R 
Dl» 
F=  D 
FRs 
FTHE 
RETU 

502 

COiMT 

ITS 

CSs 

SNs 

RXs 

Qs  R 
Vs 

ONa 
VNs 
Tls 

T2  = 

STlB 

ST2b 

F=  A 

FRs 

FTHE 

RETU 

END 

( J.K)) )  /  Dl 


OMPUTE  THE  LOCATION  CF  THE  FIXED  BOUNDARY, 

INNER  BOUNDARY, 

OUTER  BOUNDARY, 

SE  THE  EQUATION   Fs  X  *♦  lEND  +  Y  ♦*  IBND  • 

NCON 
L  .EQ.  0)EFFsNCONlN 
IBND  .GT,  DGO  TO  502 
=  COS  (THETA  (J,  K)) 
R  <J,  K) 

X  **  2  -  2,0  *  OFFSET  *  RX  *  COST  ♦  0FFSET2 
SORT  (D> 
1  -  EFF 

(RX  n,    OFFSET  *  COST)  /  Dl 
TAs  (OFFSET  ♦  RX  *  SIN  (THETA 
RN 

INUE 

IBND  •'  1 

COS  (THETA  (J,  K)) 
SIN  (THETA  (Ji  K>> 
R  (J,  K) 
X  ♦  CS  -  OFFSET 
RX  *  SN 
Q  ♦•  IT 
(V  **  IT)  * 
Q  *  QN 
V  ♦  VN 

SIGN  (1.0, 

SIGN  (1,0» 
BS  (Tl)  ♦  ABS  (T2>  -  EFF 
IBND  *  (CS  *  STl  ♦  QN  ♦  SN  * 
TAS  IBNn  ♦  HX  ♦  (-SN  *  STl  * 
RN 


NCONsQ 


PN 


Tl> 
T2) 


ST2  *  VN) 

QN  t  CS  ♦ 


ST2  •  VN) 


99 

S 

101 


MPl, 
AlO  /  IH*  ) 


SUBROUTINE  PUTOUT 
PRINT  ARRAY 
♦CALL,VARDIM 

COMMON  /SIZE/  M,  Ml, 
PRINT  99,8 
FORMAT  (//IH 
K  s  M 

PRINT  ini,  ( 
FORMAT  ((IH 
PRINT  102 
Kb  K  -  1 
IF  (  K  ,GT. 
RETURN 
102   FORMAT  (  //) 


(A,  8) 


A 


(J,  K), 

bFlS.lO 


0)  GO  TO  5 


)) 


N,  Nl,  NPl 


1,  NRl) 
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END 


SUBR 

T 

♦CAUL.LSTS 

•CALL.CBLK 

COMM 

REAL 

REAL 

REAL 

COMH 

COMM 

DIME 

EOUI 

PRIM 

FORM 

PRIN 

FORM 

CALL 

PRIN 

FORM 

IF  ( 

FORM 

FORM 

CONT 

PRIN 

FORM 

XMAX 

YMAX 

DO  9 

X 

Y 

X 

Y 

I 

I 

95       C 

( 

( 


OUTINE  PONE  (JJ) 

HE  MAIN  OUTPUT  ROUTINE, 


10 
H 


6001 

700 
701 
702 

15 


99999 


IF 
IF 

PRIM 

FORM 

CALL 

1    X 

PRIN 

99998  FORM 

PRIN 

18    FORM 

1    1 

PRIN 


ON/LA 
LAMB 
LDER 
NDER 
ON/SI 
ON/CO 
NSION 
VALEN 
T  10 
AT  (1 
T  11, 
AT  (* 
SECO 
T  600 
AT  (* 
OFFSE 
AT  (♦ 
AT  (* 
INUE 
T  15, 
AT  (1 
=  0.0 
sO.O 
5  J» 
(J) 
<J> 
(J) 
(J) 
(AB 
(AB 
ONTIN 
XMAX 
YMAX 
T  999 
AT  (1 
PLTl 
1,  J, 
T  999 
AT  (1 
T  18 
AT  (/ 
2X,  1 
T  401 


M/LAMBDA 

DA 

8,  NDERB 

T,  LDERT 

ZE/M,  Ml,  MPl,  N,  Nl,  NPl 

NST/CONST 

XI  (2),  X2  (2),  Yl  (2).  Y2  (2) 
CE  (XI,  HI),  (X2,  H2),  (Yl,  H3),  <Y2,  H4) 

HI,* *,/♦  ERRnR  SATISFIED,*) 

JJ 

ITERATIONSs*,  16) 

ND  (TIME) 
1,  TIME 

I  HAVE  BEEN  WORKING  FOR*,  rfl.S,  ♦  SEGONOS,*) 
T  ,EQ.  0,0)  CALL  OVERLAP 

FIXED  BOUNDARY*) 

FREE  BOUNDARY*) 

LAMBDA 
nx,  *LAMBDAs*,  E15,5) 


a 

3 

s 
s 

UE 
.G 

,G 
99 

Hs 
{ 
Y 

98 

HI 


N 
R  (J, 
R  (J, 
R  (J, 
R  (J, 
(X2 
(Y2 


1) 
1) 
M) 
M) 
(J)) 
(J)) 


COS 

SIN 

COS 

SIN 
,GT.  Xf'AX) 
,GT.  YMAX) 


(THETA 
(THETA 
(THETA 
(THETA 


<J, 
(J, 
(J, 
(J, 

XMAXs 
YMAX3 


D) 
D) 
M)) 

M)  ) 
ABS 
ABS 


(X2 
(Y2 


(J)) 
(J)  ) 


T,  YMAX)  YMAXs  XMAX 
T,  XMAX)  XMAXS  YMAX 


) 

IE,  T, 
1,  X2, 

) 


2, 

J, 


6,  -XMAX, 
Y2) 


XMAX,  *YMAX,  YMAX, 


*  NUM*,  12X,  IHP,  7X,  9X, 
HY,  7X,  5X,  *B  SQLAREC*) 


5HTHETA,  6X,  i?X,  IHX,  7X, 
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401   FORMAT  <*  FREE  BOUNDARY, •> 


20 
35 

402 


00 


1 

2 


1 

35  Js  2,    N 

R23  1,0  ♦  KAY2  •  R  (J,K)  **  2 
B28  C  (J,  K)  ♦  (NDERfl  {PSI,J,K)  ♦* 
(R,J#K>  *♦  2)  ♦  (NDER6  (THETA,J,K) 
♦  ZETA2 
NDERa  (PSIi  J#  K) 
TU«  LDERB  (THETAi  J,  K) 
NDER8  (THETA,  J,  K) 
LDER9  (R,  J,  K) 
NOERB  (R,  J,  K) 
RU  ♦  TV  •  RV  ♦  TU 
((PV  /  XJ)  ♦•  2)  * 
R2  /  R2 

J,  R  (J,K>,  THETA 
13,  5F20,8) 


2) 


/  (C 
2)) 


(J,K)  ♦  (NOERB 


PVa 

TVs 
RU» 
RV« 
XJa 
B28 
62  = 

PRINT  20, 
FORMAT  (IH 


(TU  ♦♦  2  ♦  C  (J,K>  ♦  RU  ♦•  2)  ♦  ZETA2 


(J#K),  XI  (J>i  Yl  (J),  82 


CONTINUE 

PRINT  402 

FORMAT  (*  FIXED  BOUNDARY,*) 

Ks  M 

DO  45  J«2,  N 

R2s  1,0  ♦  KAY2  *  R  (J,K)  *♦ 
b2=  C  (J,  K)  ♦  (NDERT  (PSI, 


2 

JiK) 


** 

1  (R#J,K)  *♦  2)  *  (NDERT  (THETA, J, K) 

2  ♦  ZETA2 
PVa  NDERT 
TUs  UDERT 
TVs  NDERT 
RUs  LDERT 
RV=  NDERT 
XJs  RU  * 
B2=  ((PV 
82=  B2  / 
PRINT  20, 

45       CONTINUE 

CAUL  CAN 

Ks  Ml 

SUMlr  SUM2S  0,0 

DO  467  J=2,  N 

R  (J*1,K)  -  R  (J-1,K) 
R  (J,K*1)  -  R  (J,K«1) 
THETA  (J*l,  K)  -  THETA 
THETA  (J,K*1)  '  THETA 


2) 


/  (C 
2)} 


(J,K)  ♦  (NDERT 


(PSI,  J,  K) 
(THETA,  J,  K) 
(THETA,  J,  K) 
(R,  Ji  K) 
(Rf  J«  K) 

TV  "  RV  *  TU 

/  XJ>  ♦♦  2)    * 

R2 

I  R  (J,K),  THETA 


(TU  *♦  2  ♦  C  (J,K>  *  RU  *♦  2)  ♦  ZETA2 


(J,K),  X2  (J),  Y?  (J),  82 


PU  = 
RVa 
TU  = 
TVs 
PU  = 
PVs 
XJ  = 
XKR: 


(J-1,K) 
(J«K-1> 


PSI  (J*1,K>  -  PSI 
PSI  (J,K*1)  -  PSI 
RU  ♦  TV  -  RV  ♦  TU 
KAY  *  R  (J,K) 


(v-l,K) 
{w,K-l) 


SIGMA2S  1,0  ♦  XKR 
Sl=  ((-(PU  *  TV  - 

*  R  (J,K)  *  TU 
S2=  (-PU  ♦  RV  ♦  PV 


PV 


TU)/XJ  ♦  X<R  ♦  Z6TA)/SIQMA2) 
*  RU)  •  RU/  (XJ  ♦  R  ( J,K) ) 
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SUMln  SUMl  ♦  SI 
SUM2a  SUM2  ♦  S2 
467      CONTINUE 

SUia  -(SUMl  ♦  SUM2  )  /  ?,0 
PRINT  469,  SUM 
469   FORMAT  (*  THE  CURRENT  INTEGRAL  IS  *,  F20,10) 
RETURN 
END 


SUB 
COM 
DIM 
♦CALL.CBU 
PR  I 


860 


861 


520 


FOR 
PRI 
FOR 

I  */ 

IQ« 

ITE 
ITE 
ITE 
ITE 
DO 


CFs 


851  FOR 

852  FOR 


ROU 

MON 

ENS 

K 

NT 

MAT 

NT 

MAT 

/) 

ITE 
R  ( 
R  ( 
R  ( 
R  ( 
890 
K  =  M 
RU  = 
PVs 
TUs 
TVs 
Ss 
Rl  = 
1. 
AF  = 
REF 
CNF 
Cl  = 
C2  = 
C3s 
C4  = 
PRI 
MAT 
PRI 
MAT 
IF 
1 

RUs 
RV  = 
TUs 


TINE  CRMANN 

/SIZF./M,  Ml,  MPl,  N,  M,  NPl 
ION  ITER  (5) 

660 


(*  C-R  EQUATIONS,  DERIVATIVES.  AND  REAL  (F)  IMAG  (F)*) 
861 

(*  SIGMA, RU  R.THETAV  SIGMA, HV  R.THETAU  RU  RV  THETAU  ThETAV 

/  4 
R  (l)s  N 
2)s  1  +  IQ  /  2 
3)s  1  ♦  lU 

4)s  ITER  (2)  ♦  ITFR  (3)  -  1 
5)3  1  ♦  2  *  10 

J-    2,  N 

-  LDFRT  (R,J»K) 

-NOERT  (R#  J.  K) 

-LDERT  (THbTA,  J,  K) 

-NDERT  (ThETA,  J,  K) 
SIGMA  (J,  K) 

R  (J,  K) 
0  /  Rl  ♦*  2 

1,0  /  S  ♦♦  2 
s  (AF  ♦(TU  ♦*  2  -  TV  •*  2)  ♦  C!"*(RU  *♦  ?  -  RV  **  2))/4,0 
B  -(AF  *  TU  *  TV  ♦  CF  *  FU  *  Rv)  /  ?,0 

S*  RU 

Rl  *  TV 

S  *  RV 

Rl  ♦  TU 
NT  B51,  J»  K,  CI,  C2,  C3,  C4,  ?U,  RV,  TU,  TV 

(♦  POINT*,  215/  Ih   ,  8F15,6  ) 
NT  852,  REF,  CNF 

(*  REAL  f-*,    F20.10,  *  IMAGINARY  F»*,  F2n,tO/) 
(K  ,NE.  M)  GO  TO  890 

LDERb(R,J,K> 
NDERP(R, J,K) 
LDERB(TH6TA, J,K) 
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TVs  NDERB(TH£TA,J,K) 

GO  TO  520 
890      CONTINUE 
RETURN 
END 


699 


SUBROUTINE  CAN 
CHECK  ANSWERS  VERSUS  THE  THEORETICAL  ANSWERS. 
CALL.CBLK 

COMMON  /SIZE/  M,  Ml,  MPl,  N,  Nl,  NPl 

IF  (IBND  ,GT,  2)RETURN 

IF  (OFFSET  .NE,  0,0)RETURN 

PRINT  100 

IF  (KAY)  699,  700,  699 

TEMPs  2ETA  /  KAY 

K  =  M 

Tl=  ZETA  ♦  KAY  ♦  (R  (1,  K)  **  2)  /  2,0 
T2=  -TEMP  *  AUOG  (R  (1,  K>) 
KbI 

Sl=  ZETA  *  KAY  ♦  (R  (1,  K)  *♦  2)  /  ?,0 
S2=  -TEMP  *  ALOG  (R  (1,  K)) 
Xa  Tl  -  T2 
Y=  SI  -  S2 
CAJs  X  -  Y 
Fl«  PSI  (1,  M)  -  T2 
F2a  PSI  (1,  1)  -  S2 
ALPhs  (Fl  •  r2>  /  CAJ 
CONSTs  (-Y  •  Fl  ♦  X  *  F2)  /  CAJ 
BETe  1,0  •  ALPH 
DO  10  Kal,  M 

EX=  ALPH  *  ZfcTA  ♦  KAY  ♦  (R  (1,  K)  ♦*  2)  /  2,0 

1  ♦  BET  ♦  (-TEMP)  *  ALOG  (R  (1,  K)) 

2  ♦  CONST 
EPS*  EX  -  PSI  (I,  K) 

PRs  ALPH  ♦  ZfcTA  *  KAY  *  R  (1,  K)  *  qET  ♦  (-TEMP)  /  R  (I,  K) 

P2s  (PR  ♦♦  2  ♦  ZETA2)  /  (1.0  ♦  KAY2  *  R  (1,K)  ♦*  2) 

PRINT  200,  K,  PSI  (1,K),  EX,  EPS,  R2 
FORMAT  (3X,  13,  F15,9,  8X,  F15.9,  8X,  E8.1,  F12.6) 

CONTINUE 
FORMAT  (/♦  EXACT  ANSWERS  VS.  COMPUTED  ANSWERS,*/ 
1    5X,  ♦NUM*,  lOX,  wEXACT*,  12X,  ♦CQviPUTED*.  lOX,  *ERRO^^*) 
PRINT  300,  ALPH,  BET 

FORMAT  (*  ALPHs*,  F10,4,  ♦  BETs*,  FlQ.4) 
PRINT  301,  CONST 
301   FORMAT  (*  CONSTANTS*,  Fl0,4) 

RETURN 
700   RX=  R  (1,1) 

CS=  CCS  (THRTA  (1,1)) 

Tls  ALOG  (RX  **  2-2,0*OFFSET*RX*CS  ♦  0rFSET2) 


200 

10 

100 


300 
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709 


RXs  R  (1,  M) 

CSs  COS  <ThfeTA  (1,M)  ) 

T2»  ALOG  (RX  ♦*  2-^,0*0FFSfcT*HX*CS 

ALPHs  (C2  -  CD  /  (T2  -  Tl) 

RfcTs  (C1*TJ>  -  G2*T1)  /  (T2  -  Tl) 
DO  709  Ksl,  M 

fX=  ALPH  ♦  AUOG(R 

♦  urrseT2)  ♦  ber 

EPS=  EX  /  PSI  (1, 

PRINT  200,  Ki  PSI 

CONTINUE 
PRINT  300,  ALPH,  HET 

RETURN 
ENU 


♦  nrFSET2) 


a,K)**  2  -2,n*TrFSET  *R(liK)*COS(TrlETA{l,K)) 


(1. 


K),  EX,  EPS 


SUBPOUTINE  FROUT  (JJ) 

PRINT  SAMPLE  VALUES  OF  F  ANT  THET&, 
♦CALL.CBLK 

COMMON/SIZE/M.Ml,  MPl,  N,  M,  NPl 

COMMON/LAM/LAMBUA 

REAL  LAMBDA 

REAL  NDERB,  NDEHT 

REAL  LDERB,  LDERT 

Js  5 
K=  3 

THETAVs  NDENB  (THETA,  J,  K) 
T=  LUERT  (THETAi  J,    K) 
Tls  NDERT  (P,  J,  K> 

RVs  (Tl  /  (R  (Ji  K>  ♦  T) )  *♦  2  -  1,0 

WRITE  OUTPUT  TAPh9Q,  712,  wJ,  P  (J,  <),  TM6TA  (J,  K),  RV.  THtTAV, 
1    LAMBDA 
712   FORMAT  (IH   ,  lb,    9X,  2F15,'i,  2E15.2,  f:20.10  ) 
K=  M 

TUVa  NDERT  (THETA,  J,  K) 
CALL  FNC  (1,  F,  FR,  FTHETA,  J,  K) 
RUVs  -F  /  FR 

WRITE  OUTPUT  TAPEb9,  712,  vJ,  R  (J,  O.     THETA  (J,  K),  RJV,  TUV 
RETURN 
END 


SUf^ROUTINE  OVERLAP 

CHECK  FOR  OVERLAP, 
CO.-IMON/SIZE/M,  Ml,  MPl, 
♦CALL,CBLK 
Ks  M 


Nl,  NPl 


125 


Kls  HI 
K2SM 
3     DO  25  Jx2,  N 

IF  (THETA  (J,  K2) 

PRINT  ion,  J,  K 

100  FORMAT  {*  THETA  OVERLAP 
7        IF  (R  (J,  K)  ,GE,  R 

PRINT  101,  J»  K 

101  FORMAT  (*  R  OVERLAP  AT 
25       CONTINUE 

IF  (K  ,EQ,  2)  GO  TO 
IF  (K  .NE.  M)  GO  TO 
Kl«l 
K2»l 
K  =  2 

QO  TO  3 
767   CONTINUE 
RETURN 
END 


,GE,  THETA  (J*l,  K2))  GO  TO  7 


AT 


♦i  215) 
KD)  JO 

215) 


TO  25 


767 
767 


SUBROUTINE  NICEPUT 

OUTPUT  THE  VALUES  OF  PSI,  R,  THETA, 
COMMON/SIZE/M,  Ml,  MPl,  N,  Ml,  NPl 
♦CALL.CBLK 

PRINT  100 

100  FORMAT  (  2(9X,  ♦PSI*,  lOX,  *R*,  16X,  *THETA*)> 
DO  15  J«l,  Mpl,  2 

K  =  M 
5    PRINT  101,  K,  PSI  (J,  K),  R  (J,  K),  THETA  (J,  K), 
1    PSI  (J*l,  K),  R  <J*1,  K),  THETA  (J*l,  K) 

101  FORMAT  (IH  ,  14,  3F15,9,  5X,  3F15,9) 
K=  K  •  1 

IF  (K)7,  7,  5 
7     PRINT  102 

102   FORMAT  <///) 
15    CONTINUE 

RETURN 

END 


REAL 
*CALL,CBLK 

SIUMA= 
RETURN 

END 


FUNCTION    SIGMA<J,K) 
=    SQRTd.O    ♦    KAY2    * 


R    (J,K)    *♦    ?) 
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REAL  FUNCTinN  LDERT  (A,  ^,     K) 
ALL  SUBROUTINES  CALLINLi  TUS  FUNCTION  MUST  DECLARE  IT  TO  dfc  HEAL.,, 
♦CALL.VARDIM 

COMMON  /MESH/  H2 

COMMON  /SI2F/  M,  Ml,  MPl,  N,  Nl,  NPl 
ROUTINE  TO  COMPUTE  LINE  DERIV  ON  TOP, 

LDfcRTs  {-'A  <J  ♦  li  M)  ♦  A  (J  ^  1,  M))  /  H2 

RETURN 

END 


REAL  FUNCTION  NDtHT  (A,  ^,    K) 
ROUTINE  TO  COMPUTE  INNER  NCRh-AL  TERIV  ON  TOP, 
•CALL.VARDIM 

C0MM0N/MESH/H2 

COMMON  /SIZE/  Ml  Ml,  MPi,  N,  Nl,  NPl 
ALL  SUBROUTINES  CALLING  ThIS  FUNCTION  MUST  DECLARE  IT 
NDERTs  (-3,0  *  A  (J,  M)  ■»  4.0  *  A  (J,  M  -  1)  -  A  (J, 
RETURN 
END 


TO  tJE  REAL,  ,, 
M  -  2>)  /  M2 


REAL  FUNCTION  LD6R8  (A,  ^,    K) 
n      COMPUTE  LINE  UERIV  IN  THE  f^CTTOM, 

C        ALL  SUBROUTINES  CALLING  T(-IS  FUNCTION  MUST  OECLARE  IT  TO  dE  REAL 
*CALL»VARDIH 

C0MM0N/MESH/H2 

LDERBs  (A  (J  ♦  1,  1)  -  A  (v  »  1,  1))  /  H2 

RETURN 

END 


REAL  FUNCTION  NOERB  (A,  ^,  K) 
C   COMPUTE  INNER  NORMAL  UERIV  ON  BOTTOM, 
♦CALL.VARDIM 

ALL  SUBROUTINES  CALLING  TI-IS  FUNCTION  MUST  DECLARE  IT  TO  BE  REAL. 
C0MM0N/MESH/H2 

ND6P8=  (-3,0  *  A  (J,l>  ♦  4,n  *  A  (J, 2)  "  A  (J,3>)  /  H? 
RETURN 
END 
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SUB 
ITERA 
♦CALL.CBL 
EOU 
COM 
COM 
COM 
COM 
COM 
EMA 
Jls 
IPs 
DO  SE 
J  = 
DO 


IF 

Ds 


/  4,0 


31   CON 

MOVE 
115 
105 

HERAT 
DO 


IF 
Fs 


ROUTINE  PSINT 

TE  PSI  ON  THE  INTERIOR, 

K 

IVALENCE   (Z  (l)i  EMAX). 

MON/REL/WPP,i^RP,WTP 

HON  /SIZE/M,  Ml,  MPl,  N, 

MON  /MESH/  U2 

HON  /CNT/  IP.  ITi  IR 

HON  /TIME/  RP,  RR,  RT, 

X=  0,0 

Kla  0 

IP  *  1 
COND  COLUMN, 
2 

105  Ks  2,    Ml 
DEL=        (PSI  (J  ♦  1, 

♦  PSI  (J,  K  -  1)  ) 
(KAY  ,EQ.  0,0)  GO  TO  31 
CALCD  (J,  K> 

CALCF  (J,  K) 

TU=  THETA  (J*1,K)  -  ThETA  (J-1,K) 

TVs  THETA  (JiK*l)  -  TkETA(J,K-l> 

T=  1,0  /  SORT  (C  (J»  K)) 

r<U=  R  (J*1,K)  -  R(J-1,K) 

RV=  R  (J,K*1)  -  R  (J,K<.1) 

PSIUs  PSI  (J*1.K)  -  PSI  (J-l.K) 

PSIVa  PSI(J,K+1)  -  PSI  (J,K-1) 

DELS  DEL  ♦  (GR  (J,  K)  /  (2,0  *  C 

♦  D)  *  T  *  (TV*PSIU  -  TL*PSIV)  / 

DEL*  DEL  ♦  F  ♦  T  *  (RL  *TV  -  RV 

IF  (ABS  (DEL)  .Lfc.  ABS 

EMAXx  DEL 

Jl=  J 

Kls  K 

TINUE 

PSI  (J,  K)a  PSI  (J.  K) 

SECOND  COLUMN, 
PSI  (NPl,  K)x  PSI  (2, 
CONTINUE 
E  REST, 
125  Js  3, 
no  125  Ks 
DELS 

♦  PSI 
(KAY  ,EQ, 
CALCD  (J, 
CALCF  (J, 
TUs  THETA 
TVs  THETA 
T=  1,0  / 


(Z  (2),  Jl),  (Z  (3).  Kl) 
M,  NPl 


RFR,  RFX 


K)  ♦  PSI  (J  "    1»  K) 
(J,  K) 


♦  PSI  (Jl  K  ♦  1) 


PSI 


(EMAX)) 


(J,  K)) 

16.0 
♦  TU)  /  16,0 
GO  TO  31 


DEL 


/<PP 


K) 


N 

2,    Ml 
(PSI  (J  ♦  1, 
(J,  K  -  1)  ) 
0.0)  QH  TO  33 
K> 
K) 

(J*liK)  -  THETA  (J-1,K) 
(J.K*1)  -  TI-ETA(J,K-1) 
SORT  (C  (J»  K)) 


K)  ♦  PSI  (J  -  1,  K) 
/  4,n  -  P^I  (J.  K) 


♦  PSI  (J.  K  ♦  I) 


RUs  R 
RV=  R 


(J*1.K)  - 
(J,K*1)  - 


R( J-l.K) 
R  (J.K-l) 
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33 

125 

MOVE 
DO 
135 


PSIUs  PS!  {J*1,K)  -  PSI  (J-1,K) 

PSIV=  PSI(J#K*l)  -  PSI  (J,K-1) 

DELS  DEL  ♦  (CR  (J,  K)  /  (2.0  *  C  (J,  K)) 

♦  D)  *  T  *  (T\/*PSIU  -  TL*PSIV)  /  16.0 

PELS  DRL  ♦  F  *  T  *  (RL  *Tv  -  «V  ♦  TU)  /  16,0 

IF  (ABS  (DEL)  ,LE.  ABS  (EMAX))  GO  TO  33 

EMAXs  DEL 

Jl=  J 

Kl=  K 

CONTINUE 

PSI  (J,  K)5  PSI  (J,  K)  ♦  DEL  *  WPP 
CONTINUE 
INTO  THE  FIRST  COLUMN, 
Ks  1.  M 
(1.  K)=  PSI  (N,  K) 


135 
PSI 

RETURN 
ENU 


SU9R0UTINE  THE  I  NT 
ITERATE  THETA  IN  THb  iNTETIOfi, 
♦CALL.CBLK 

EQUIVALENCE  (Z  (7),  £MAX),  (Z  (3),  Jl),  (7  (9)«  Kl) 
COMMON  /CNT/  IPi  IT.  IR 
COMMON  /SIZE/M,  Ml,  MPl,  N,  M,  i>jPl 
COMMON  /TIME/  RP,  RR,  RT,  PFR,  RFX 
COMMON/REL/WPP,WRP,WTP 
EMAXs  0,0 
Jls  Kl=  0 
IT=  IT  ♦  1 
DO  SECOND  COLUMN 
Js  2 

UO  ?5  K=NA,  M2NA,  NA 
KP=  J  ♦  K 

DEL*  (TH  (KP*1)  ♦  TH  (KF-1)  ♦  TH  (KP+NA) 
1    -  TH  (KP) 

IF  (ABS  (DEL)  ,LE,  ABS  (EMAX))  QO    TO  31 

EMAXs  DEL 
Kl=  KP 
31   CONTINUE 


*  TH  (KP-NA))  /  4iO 


TH  (KP) 
KP)=  TH 


TH  (KP)s 
TH  (Nl  * 
?5       CONTINUE 
C       MOVE  SECOND  COLUMN, 
C       CONTINUE  ITERATIO'M 
UO  45  J«  3,  N 
DO  45  Kb.NA,  M2NA, 
KPs  J  ♦  K 
DELS  (TH  (KP*1> 


♦  DEL 
(KP) 


'  WTP 
PICCN 


NA 
*    TH  (KF-1)  ♦  TH  (KP*NA)  ♦  TH  (KP-NA))  /  4,0 


129 


33 


-  TH  (KP) 

If    (ABS  (DEL)  .LE,  ABS  (EMAX))  GO  TO  33 
EMAXa  DEL 
Kl3  KP 

C0NTINU6 


TH  (KP>r  TH  (KP)  ♦  DEL  *  WTF 
45       CONTINUE 

MOVE  INTO  THE  FIRST  COLUMN, 

MlNAls  MINA  ♦  1 

DO  55  K«l,  HlNAl,  NA 
55       TH  (K)3  TH  (Nl  ♦  K)  ♦  PICON 

Jls  MOD  (Kl,  NA) 

Kls  Kl  /  NA  ♦  1 

RETURN 

END 


SUBROUTINE  PINT 
♦CALL.CBLK 

EQUIVALENCE  (Z  (4),  EMAX), 
CCiMON  /SIZE/M,  Ml,  MPl,  N, 
COMMON/REL/WPP,WRP,WTP 
COMMON  /CNT/  IP,  ITi  IR 
COMMON  /TIME/  RP,  RR,  RT, 
ITERATE  R  ON  THE  INTERIOR, 
IRa  IR  ♦  1 
EMAXs  0.0 
Jls  0 
Kl*  0 
ITERATE  THE  SECOND  COLUMN, 
Js  2 
DO  65  Ks  2,  Ml 

C32S  32. n  *  C  (J,  K) 
f'Us  R  (J*1,K)  -  H  (J-1,K) 
RVs  R  {J,K*1>  -  R  (J,K-1) 
DELS  (R  (J  ♦  1,  K)  ♦  R  (J  - 

1  )  /  4.0  -  R  (J,  K) 

2  +  CR  (J,K)  ♦  (HU**2 
IF  (  ABS  (DEL)  .LE. 
EMAXs  DEL 

Jls  J 
Kls  K 
101   CONTINUE 


i 


(Z  (3),  Jl), 
M,  MPl 


RFR,  RFX 


(Z  (6),  Kl) 


i 


1,  K)  ♦  R  (J,  K  *  I)  ♦  R  (J,  K  -  1) 


♦  RV  **  2)/:32 

A6S  (EMAX))  GO 


TO  101 


R  (J,  K)s  R  (J,  K)  ♦ 
C  (J,  K)«  CALCC  (J,  K) 
MOVE  SECOND  COLUMN, 

R  (NPl,  K)e  R  (2,  K) 
C  (NPl,  K)c  C  (2,  K) 


EEL  *  ^RP 
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65       CONTINUE 

DO  THt  REST  CF  THb  ITfcRATlCN, 
UO  65  Ja  3,  N 

DO  B5  K=  2i  Ml 

C32S  32. U  *  C  (J,  K) 

PU=  R  (J*1,K)  -  H  (J-1,K) 

RV=  H  ( J,K+1)  -  K  ( J.K»1 ) 

DfeL*  <R  (J  ♦  1.  K)  ♦  p  (J 

1  )/4.0-t^(.J#K) 

2  *  CR  (J.K)  *  (RU**2  ♦ 

IF  (ABS  (UEU)  ,L£. 
EMAXs  DhL 
Jle  J 
Kl«  K 
103   CONTINUE 


-  1,  K)  ♦  R  (J,  K  ♦  1)  ♦  R  (J,  K 


RV  **  2)/:32 

ABS  (EMAX))  GO  TO  11)3 


1) 


R  (J,  K)a 
C  (J,  K)= 
85  CONTINUE 

MOVE  INTO  THE  FIRST 
00  95  Ks  1,  M 

R  (1,  K>B  R  <N, 
95       C  (li  K)s  C  (N, 
RETURN 
END 


R    (J/ 
CALCC 


K) 


•    DEL 
K) 


W^P 


COLUMN, 

K) 
K) 


MPl,     N,     M,     IxlPl 


ENFR 


SUapOUTINE  FREBNO 
ITERATE  FREE  BOUNDARY, 
*CALL,LSTS 
•CALL.CBLK 

EQUIVALENCE  (Z  (14),  EMaX),  (2  (1.'5),  Jl),  (Z  (1&).  Fl),  (Z  (l7> 
I    Jie) 

COMMON  /SIZE/M,  Ml, 

COMMON  /HND/  IFX,  IFR 

COMMON  /b/BTFX,  «F,  8TFR, 

COMMON/LAM/LAMBUA 

COMMON  /M6SH/DNM 

REAL  LAMBDA 

DIMENSION  RUS  (?)•  RVS  (2),  THETAUS 

EQUIVALENCE  (RUb,  HI),  (RVS,  H2 ) 

EQUIVALENCE  (THETAUS,  M3),  (THETAVS 

EQUIVALENCE  (PSIV,  H5) 

REAL  LDERT,  NDERT,  LDERB,  NUERE 

COflMON/FXCON/    JJ,     KK 

DIMENSION  ST  (7,  1«) 

COMMON  /INT/TOPR,  TUPT,  EOTR,  EOTT 

DIMENSION  TOPR  (140),  TOFT  (140),  BOTR 

SIQMA2  (J,K)s  1,0  *  KAY2  *  R  («,  K)  »* 

EMAXs  0,0 


(2),  THETAVS  (2),  pSIV  (2> 


H4> 


(140), 
2 


B  0  T  T  (14  0) 
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25 


0)  GO  TO  3 


N 

LDErtB 
i  NDERH 

(  J)8 


LDERB 

NDERB 


J,  K) 
J,  K) 
(TI-ETA, 
(ThETA 


J.  K) 
J.  K) 


El»  0,0 
Jlx  J28  0 

IF  (IKON  ,LE, 
CALL  CONFORM 
GO  TO  47 
CONTINUE 
COMPUTE  LAMBDA 
SUMr  0,0 
SUMPS  0.0 
K=  1 

DO  25  J«  2, 
PUS  (J)s 
RVS  (J)= 
THETAUS 
THETAVS 
PSIV  (J)B  NDtRb  (PSI,  J.  K) 
T3s  C  (J,  K>  ♦  RUS  (J)  **  2  ♦  THETAUS 
T7Ps  (C  (J,K)  *  PSIV  (J)  *♦  2)    /    T3 
T7=  (T7P  ♦  ZETA2)  /  SIQh'Aa  (J,K) 
ERATUs  2FTA2  /  T7 
SUMP*  SUMP  ♦  T7 
SUMS  SUM  +  (PSIV 
CONTINUE 
LAMBDAS  SUM  /  Nl 
TLAM=  SUMP  /  Nl 
JSs  MIN0(N,16) 
UO  45  J=  2,  N 

T=  C  (J,  K)  *  RUS  (J) 

C  (J,  K)  ♦  RUS  (J)  * 

/  T 

PSIV  (J)  **  2 

C  (J,  K)  *  RVS 

T3  *  ZETA2  /  C 

LAMBDA  ♦ 


(J)  **  2 


(J)  *♦  2)  /  T3  ♦  ZETA2  /  C  (J,  K) 


** 


Tl  = 
T4s  Tl 
T2  = 
T3  = 
T8  = 
T5s 
T2s 
T6s 
ERATV 
TlOOs 


2  ♦ 
RVS 


THETAUS  (J)  **  '^. 
(J)  ♦  THETAUS  (J) 


•  THETAVS  (J) 


C  (J*  K) 

(T2/T3  ♦ 

ZETA2 

TLAM  / 


(J) 
(Ji 
T3  /  (T2 
♦  PSIV  (v 


*♦  2  ♦  THETAVS  (J)  *♦  2 

K) 
♦  T8)  • 
>  **  2 


1.0 


ZETA2) 
/  T6 

T6  -  1.0 
16) 
T5 
ST  (2,  JT)s  iJs  TLAM  - 
ST  (3,  JT)=  U  /  T6 
ST  (4,  JT)s  CJ  /  TLAM 
ST  (5,  JT)s  Q    /    JS 
RATS  SORT  (RVS  (J) 
ST  (6,  JT)=  ST  (3, 
ST  (7,  JT)s  ST  (4, 
IXs  IA8S  (ISWT) 
T5s  ST  ( IX,  JT) 
IF  (ABS  (T4)  .LE, 
EMAXs  T4 
Jls  J 


jTs  MINO  (J. 

ST  (1,  JT)= 
(2,  JT)s 
(3,  JT)= 
(4,  JT)s 
(5,  JT)s 


/  SIfiMA2  (J,K) 


T6 


JT) 
JT) 


2  f  THETAVS 
/  RAT 
/  RAT 


(J)  **  2) 


ABS  (EMAX))  QO  TO  31 


132 


31 


El 


33 


CUNTINUE 
IF  (  ARS 
i  T5 
J2=  J 
CONTINUE 


(T5)  ,l&.  AHS  (fil))  GO  n  33 


T4WS 
T5^^a 

BOTH 
P  (J, 
C  (J, 

DELS 
BOTT 


biTFR 
tiNFH 
T4W 
(J)  = 
K)s 

T4W 
(J)s  DEL 


♦  T4 

*  TiJ 

*  RUS)  (J)  ♦  T5^    *  RVS  (J) 
DEL 

R  (J.  K)  ♦  CEL 
CALCC  (J,  K) 

♦  THfcTAUS  (J)  ♦  T5U  *  TH5TAVS  (J) 


THETA  (J, 
45       CONTINUE 
LAMBDAS  TLAM 
MOVE 
47    R  (l»  1)=  R  (N, 
THETA  (l,l)s 
C  (1,  l>s  C  (N, 
R  (NP1» 

THETA  (NPl, 
C  (NPl,  l)s  c 
RETURN 
END 


K)»  THETA  (.^,K)  ♦  DEL 


X) 

THETA 
1) 
1)=  R  (2,  1) 


(N,l)  *  FICON 


1)=  THETA  (2,  1) 
(2,  1) 


P120H 


SUBROUTINE  FIXHND 
ITERATE  FIXED  BOUNDARY 
♦CALL.CBLK 

EQUIVALENCE  (Z  (10),  EMAX),  (Z  (11),  Jl),  (Z  (12).  Fl),  U    (l3), 
1    J2) 
EQUIVALENCE  (Tl,  OMEGA) 
COMMON  /SIZE/M,  Ml,  MPl,  N,  M,  gPl 
COMMON  /B/BTFX,  BF,  HTFR,  BNFR 
REAL  LDERT,  NDEHT,  LDERB,  NDEHF 
COMMON  /INT/TOPR,  TORT,  EOTR,  BOTT 

DIMENSION  TOPR  (140),  ToFT  (140),  HOTR  (140),  niOTT  (140) 
DATA  PI/3, 14159265/ 
L  =  l 

EMAX=  0.0 
El*  0.0 
Jl*  J2=  0 
K=  M 

IFXe  IFX  ♦  1 
J=  ? 

BU=  LDFRT  (R,  J,  K) 

RVs  NUERT  (R,  J,  K) 

THETAUa  LDERT  (THETA,  J,  K) 


133 


131 


133 


C 
I 
El« 

J 
C 
T 
T5Ws 
D 
T 

C 

D 

T 

T 

R  (N 


(N 


31 


33 


THei 

DO  2 

P 
R 
T 
T 
C 
T 
T 
T 


C 
I 

tla 
J 
C 
T 

T5W  = 


HETAV 
ALL  F 
9=  C 
4=  (C 
5=  F 
IF  ( 
EMAX 
Jl« 
ONTIN 
F  <  A 
T5 

2=  J 
ONTIN 
4WS  B 
T5  * 
EL«  T 
OPR  ( 
(J. 
(J. 
EL=  T 
OPT  ( 
HETA 
PI*  K 
Pl#  K 
A  (NP 
05  Ja 
U=  LD 
V=  ND 
HETAU 
HETAV 
ALL  F 
9s  C 
4s  (C 
5s  F 
IF  ( 
EMAX 
Jl  = 
ONTIN 
F  (  A 
T5 

2s  J 
ONTIN 
4W=  B 
T5  * 
EL=  T 


s  NDERT  ( 
NC  <L.  F, 
(J,  K)  * 

(J,  K)  * 
/  (FR  •* 
ABS  (T4) 

s  T4 

J 

UE 

BS  (T5)  .LE,  ABS  (ED)  GO  TO  l33 


THETA,  J,    K) 

FR,  FTHETA,  J,  K) 
RU  *•  2  ♦  THETAU  **  2 

RU  *  RV  -♦•  ThETAU  ♦  THETAV)  /  T9 
2  ♦  FThETA  ♦♦2) 
,LE,  ABS  (EMAX))  GD  TO  131 


*  T4 


UE 
TFX 

BF 
4W  ♦  RU  - 
J>B  DEL 
K)8  R(J, 
K)a  CALCC 
4W  ♦  THET 
J)«  DEL 
(J,  K)=  T 
)s  R{2.  K 
)s  C  (2i 
1,  K)=  TH 

3,  N 
ERT  (R, 
ERT  (R, 
s  LDERT 
s  MDERT 
NC  (L, 


J 
J 
( 
( 
F, 


T5W  *  FR 

K)  ♦  DEL 
(J#  K) 

AU  •>    T5W  *  FTHETA 

HETA  (w»  K)  ♦  DEL 

> 

K) 

ETA  (2,  K)  •  PICON 

,     K) 
,    K) 

THETA,  J,  K) 
THETA.  J,  K) 

FR,  FTHETA,  J,  K) 
RU  **  2  ♦  THETAU  **  2 

HU  *  PV  ♦  THETAU  *  THETAV)  /  T9 
2  ♦  FTKETA  ♦*2) 
,LE,  AES  (EMAX) )  GO  TO  31 


(J,  K)  * 

(J,  K)  * 
/  (FR  ** 

ABS  (T4) 

s  T4 

J 

UE 

BS  (T5)  ,LE,  ABS  (ED)  c90  TO  33 


UE 
TFX 
BF 


T4 


RU 


T5W  *  FR 


TUPR  (J)s  DEL 
R  (J,  K)s  R(J,  K)  + 
C  (J,  K)B  CALCC  (J, 
DELS  T4W  ♦  THETAU  - 


DEL 
K) 

T5W 


*  FTHETA 


TOPT  (J)=  DEL 


134 


20  5 


C 

c 
c 


520 


521 


570 


571 


700 

I 
701 

703 
522 


305 
PU 


T 

C 

MOVt 

R  (1 

C  (1 

THfeT 

N 


IF  { 
IF(I 
J=  < 
Jls 
JP1« 
K=  h 
GO  T 
CONT 
IF  ( 
K=  M 
«0  T 
CONT 
IF  ( 
Ka  ( 
CONT 
J=  N 
Jl« 
JP1> 
CONT 
IF(J 

U18 

IF  ( 
Q2  = 
IF  ( 
PRIN 
FORM 
RfcTU 

OOK    T 
Tl  = 
CALL 
RETIJ 

LOOK 
Tls 
CALL 
RETU 
CONT 
INTE 
SUHs 
DO  3 
S 

T  IN 
Ts  S 
SUMb 


HETA  (J,  K)=  THhTA  (v.  K)  *  DEL 
ONTINUfe 

,  K)s  R  (N,  K) 
,  K)=  C  (N,  K> 
A  (1,  K>a  THfcTA  (N,  K)  *  PICON 

ormalizing  schemes, 
ipoinTbo  none 

OTHERWISE  THETAsU  aT  ONE  POINT. 
IPOINT  ,EQ,  0)RhTURM 
POINT  .NE,  10)GO  TO  520 
NPl  *  1)  /  2 
J  -  1 

J  *  1 

0  571 

INUE 

IPUINT  ,NE,  3)  UO  TO  521 

0  570 
INUE 

IPOINT  .NE.  4)  fiO  TO  522 

M*l)  /  2 

INUE 

Nl 

NPl 
INUE 

POINT  .NF,  iJ)Ksl 
THETA  (Jl,  K)  ♦  THETA  (J,  K) 

01  .LE.  0,0)  UO  TO  701 
THETA  (J,  K)  *  THETA  (wPl,  K) 
(J2  ,LE,  0,0)  (iO  TU  703 

T  700 

AT  (*  I  CANNOT  FIND  THETA  EQUAL  ;?cRO,*) 

RN 

0  LEFT 

THETA  (J,K)  /  (THETA  (w,K)  -  TMETA  (Jl,K)) 

PSHIFT 
RN 

TO  RIGHT 
THETA  (J,K)  /  (THETA  (v#K)  -  THETA  (JPl,  K)) 

LSHIFT 
RN 

INUE 
URATION  STARTS  HERE. 

0,0 
05  Js  2i  N 

UMs  SUM  ♦  THETA  (J,  K) 
PI  FOR  HOUSEKEEPING, 
UH  ♦  PI 

T  /  Nl 
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405 


DEL"  SUM  -  PI 

ir  (IPOINT  .EO,  2)DeL=  THETA  (N,  M) 

DO  405  Ja  1,  NPl 

DO  405  Kal,  M 

THETA  (J,  K)»  THETA  (J,    K)  -  DEL 
RETURN 
END 


SUBR 

C  D 

♦CALL.CBLK 

EQUI 

1    J 
COMM 

REAL 
COMM 
Els 

Jlaj 
L  =  0 
K  =  l 
DO  1 
P 
RVs 

THET 
THET 
CAUL 
T 
T 
T 
T 
T 
I 


131 


133 


150 


OUTINE  CONFORM 

0  A  CONFORMAL  MAPPING, 

VALENCE  (Z  {10)»  EMAX),  (Z  (11),  Jl),  (Z  (12),  Fl),  (Z  (l3), 
2) 

ON/B/BTFX,  HF,  BTFR,  BNFR 

LDERT,  NDEKT,  LDERB,  NPERB 
0N/SIZE/M,M1,MP1,  N,  Nl,  NPl 
E2«0,0 
2  =  0 


50  J 
U=  L 
NDER 
AUs 
AVs 

FNC 
9=  C 


4s 
5s 
Is 
2s 

F 


C 
F 
T 
F 
(A 


I 


T 

T 

D 

R 

C 

D 

T 

EMAX 

RETU 

END 


Els 

Jl  = 

F  (A 

E2s 

J2s 
CON 

IWa 

2Wa 

EL  = 
(J, 
(Jl 

bL  = 

HETA 

s  £2 

HN 


=  2, 

DER 
B  ( 
LDE 
NDE 

(L 

(J 

(J 
R  ♦ 
4  / 

/ 
BS 

Tl 
J 
BS 

T2 

J 
TIN 
BTF 
8NF 
TIW 
K)  = 
K)s 

TIW 
(J 


Ri  J,  K) 
Jl  K) 

(THETA,  J,  K) 
(THETA,  J,  K> 
FR,  FTHFTA, 
)  *  RU  **  2  ♦ 
)  *  RU  *  RV  ♦ 
♦  FTHETA  **  : 


**  2 

*  THETAV 


N 
B    ( 
R, 
RB 
RB 

,  F,  FR,  FTHFTA,  J,  K) 
,  K)  *  RU  **  2  ♦  THETAU 
,  K)  *  RU  *  RV  ♦  THETAU 
♦  2 

T9 
T5 
(El)  ,UT,  ABS  (Tl)  )  GO  TO  131 


(E2)  ,GT.  AB<?  (T2))  GO  TO  133 


IJE 

R  *  Tl 

R  ♦  T2 

*  HU  -  T2W  ♦  FR 
R  (J,  K)  ♦  DEL 
CALCC  (J,K) 

*  THETAU  "  T2W  ♦  FTHPTA 
,K)s  THETA  (J,K)  ♦  TEL 
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SUBROUTINE  . 
NORMALIZE 


*CALL»CbLK 

COMMOK 


I  SHIFT 

THfc  SOLUTION 


Q  = 

DO 
DO 


/SIZE/  M,  Ml,  MPl,  N,  Nl,  MPl 
OMEGA 
1,0  -  OHPGA 


* 


5  jsl,  N 

5  Ksl,  M 

PSI  (J,K)= 

R  (J,K)=  D 

THETA  (J,K)s  n 

C  (J.K)s  CALCC 

CONTINUE 

DO  7  K=l, 

PSI  (NPl, 

P  (NPl,  k)a  R 

THETA  (NPl, 

C  (NPl,  K)s 

CONTINUE 

RETURN 

END 


PSI  (J,K) 
(J,K)  ♦  G 


*  THETA 
(  J,K) 


(J,K) 


*PSI 
(  J*l, 

*  G  » 


( J+l.K) 
K) 

THETA  (Jfl.K) 


M 


PSI  (2,K) 

(2,  K) 
K)s  THETA  (Z,  K)  -  PI:ON 
C  (2,  K) 


SUBROUTINE  RSHIFT 

NORMALIZE  THE  SOLU 
COMHON  /SIZE/  M,  Ml, 
*CALL,CPLK 

Qs  CMEGA 

1,0  - 

NPl 

5  K5l,  M 

PSI  (J,K) 


TION, 
MPl,  N, 


Nl,  NPl 


D! 

UO 


OMEGA 


5  K5l,  M 

PSI  (J,K)s  D  *  PSI  (J,K)  ♦  G  *  PSI  (J-: 

P  (J,K)s  D  *  H  (J,K>  ♦  G  *  P  (J-1,K) 

THETA  (J,K)=  D  *  THETA  (J,K)  *  G  »  THE' 

r  (J,K)s  CALUC  (J,K) 

CONTINUE 


i,K) 

TA  (J-1,K) 


CONTINUE 
'  -  1 

.GT.  1) 

M 


J=  J  -  1 

IF  (J  .GT.  1)  GO  TO  3 

00  7  K=l,  M 

PSI  (1,  K)s  PSI  (N,K) 

P  (l,K)s  H  (N,K) 

THETA  (1,  K)s  THETA  (S,  K)  *    PICOM 

C  (1,K)=  C  (i^,K) 

CONTINUE 
RETURN 
END 
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EOR 


SSWITCH 

SPARl  $ 

SPARl 

$PAR2 

$PAR3 

$PAR3 

$PAR3 

$PAR3 


S 

MlBlO*  Nls34S 

$ 

1=7, 

I»6$ 

I»7, 

i»n$ 


11=10,  ERR3l,nE-4$ 

11=1,  I2al,  I33l,  ERR3l,0E-8,  iKONs-lJ 


EOR 


IF 


15 


101 


,0 
.0 


Js  1,  20 
TAPE65,  X(J) 


ADJUST 
ADJUST 
ADJUST 
ADJUST 
ADJUST 
ADJUST 


Y2(J>,  Y3(J>,  Y4(J),  Y5<J),  Y6(J) 


PROGRAM  RESID  (OUTPUT,  TAPE65) 

DIMENSION  X  (500),  YKSOO),  Y2  (  500  > ,  y3  (  500  ) ,  Y4  (  500  > .  Y5  (  500  ) 

DIMENSION  Y6  (500) 

NTs  65 

B=  0 

S»  1 

DO  0 

READ 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

IF  (ENDFILE 

CONTINUE 
NO  ENDFILE  KEEPREADING. 
DO  15  J=  21,  500 
READ  TAPE65,  X(J), 
ADJUST 


Yl 

Y2 

Y3 

(Y4 

<Y5 

(Y6 


»  YKJ), 

(J)) 

(J)) 

(J)> 

(J)) 

(J>) 

(J)) 


NT)  25,  « 


CALL 

CALL  ADJUST 
CALL  ADJUST 
CALL  ADJUST 
CALL  ADJUST 
CALL  ADJUST 

IF  (ENDFILE 
CONTINUE 
J=  500 
PRINT  101 
FORMAT  (♦ 
QO  TO  25 


YKJ), 
(YKJ)) 
(Y2(J)) 
(Y3(J)) 
(Y4  (J) ) 
(Y5  (J)) 
(Y6  (J>) 
NT)  25, 


Y2(J),  Y3{J).  Y4(J),  Y5(J),  Y6(J) 


15 


TOO  MANY  PLOTS. ,.*) 
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Vil-    HAVE  WEACHSn  AN  bND  OF  FILE, 
25    J=  J  -  1 

IF<J  .fcO,  0)  GO  TO  27 
Z=  X  (J) 
PRIKT  102 
102   FOKMAT  (IHl  ) 
S=  0,0 
bs  ?0,0 
PRINT  99999 
99999  FOWMAT  (1H=  ) 

CALL  PLTl  ( IE,  T,  1,  5,  S 
1    Yl,  J,  X) 

CALL  PLTl  (  IE,  T,  1,  b,    S 
1    Y2,  J,  )() 

CALL  PLTl  ( IE,  T,  1,  b,    S 
1    Y3,  J,  X) 

CALL  PLTl  (  IE,  T,  1,  5,    S 
1    Y4,  J,  X) 

CALL  PLTl  (IE,  T,  1,  5,  S 
1    Y5,  J,  X) 

CALL  PLTl  (IE,  T,  1,  b,  S 
1    Y6,  J,  X) 
PRINT  201,  IE 
201   FORf^AT  (♦  ERRORS*,  14) 
27    CALL  EXIT 
END 


B,  o.n,  z. 

^,  0,0,  z, 

^^  0.0,  Z. 

H,  0,0,  z. 

H ,  0,0,  Z . 

^',  0,0,  Z, 


EOF 


SUBROUTINE  ADJUST  (e) 

IF  (ABS  (E)  ,LT,  l,')E-in)  6si,0F-10 

T=  -SIUN  (1,0,  b)  •  ALOGIO  (ABS  (E)) 

IF  (T  ,LT,  0.0)Ts  T  ♦  20,0 

E=  T 

RETURN 
END 
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